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ABSTRACT 

We investigate galaxy bias in the framework of the "Halo Occupation Distribution" 
(HOD), which defines the bias of a population of galaxies by the conditional probability 
P(N\M) that a dark matter halo of virial mass M contains N galaxies, together 
with prescriptions that specify the relative spatial and velocity distributions of 
galaxies and dark matter within halos. By populating the halos of a cosmological 
N-body simulation using a variety of HOD models, we examine the sensitivity of 
different galaxy clustering statistics to properties of the HOD. The galaxy correlation 
function responds to different aspects of P{N\M) on different scales. Obtaining 
the observed power-law form of (,g{r) requires rather specific combinations of HOD 
parameters, implying a strong constraint on the physics of galaxy formation; the 
success of numerical and semi-analytic models in reproducing this form is entirely 
non-trivial. Other clustering statistics such as the galaxy-mass correlation function, 
the bispectrum, the void probability function, the pairwise velocity dispersion, and the 
group multiplicity function are sensitive to different combinations of HOD parameters 
and thus provide complementary information about galaxy bias. We outline a strategy 
for determining the HOD empirically from redshift survey data. This method starts 
from an assumed cosmological model, but we argue that cosmological and HOD 
parameters will have non-degenerate effects on galaxy clustering, so that a substantially 
incorrect cosmological model will not reproduce the observations for any choice of 
HOD. Empirical determinations of the HOD as a function of galaxy type from the 
2dF and SDSS redshift surveys will provide a detailed target for theories of galaxy 
formation, insight into the origin of galaxy properties, and sharper tests of cosmological 
models. 



Subject headings: cosmology: theory, galaxies: formation, large-scale structure of 
universe 



1. Introduction 



The relation between the galaxy and dark matter distributions depends on the physics of 
galaxy formation, and it is expected that galaxies are, at least to some degree, biased tracers of 
the mass distribution. This expectation, which is supported by observational evidence that galaxy 



clustering varies with luminosity, morphology, and color (Guzzo et al. 1997; Norberg et al. 2001 
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Zehavi et al. 2001; and references therein), complicates efforts to test cosmological models against 
observed galaxy clustering. However, the presence of bias also implies that galaxy clustering can 
be used to constrain the physics of galaxy formation, especially as independent observations define 



the background cosmology with increasing precision (e.g., Wang, Tegmark, & Zaldarriaga 2001). 



The 2dF and SDSS galaxy redshift surveys ( |Colless et al. 2001 ; York et al. 2000 ), which can 



measure the clustering of different galaxy types with unprecedented detail, are now bringing this 
goal within reach. Achieving it requires a language for describing bias that is powerful enough 
to capture the information in these measurements and thereby connect observations of galaxy 
clustering to the physics of galaxy formation. 

In this paper, we examine the influence of bias on galaxy clustering statistics, using the 
framework of the "Halo Occupation Distribution" (HOD). This approach describes bias at the level 
of "virialized" dark matter halos, structures of typical overdensity p/p ~ 200, which are expected 
to be in approximate dynamical equilibrium. Gas dynamics, radiative cooling, and star formation 
can strongly influence the distribution of galaxies within such halos (e.g., the numbers, masses, 
and locations of galaxies), but the masses and spatial distribution of halos themselves should be 
determined mainly by gravitational dynamics of the dark matter. In the HOD framework, the bias 
of any particular class of galaxies is fully defined by the probability distribution P{N\M) that a 
halo of virial mass M contains galaxies, along with the relations between the galaxy and dark 
matter spatial and velocity distributions within halos. While the history of the galaxy population 
is necessarily entwined with the background cosmology, the HOD description suggests a useful 
conceptual division between the "cosmological model" and the "theory of galaxy formation" in 
predictions of galaxy clustering: the cosmological model determines the properties of the halo 
distribution, and the theory of galaxy formation specifies how those halos are populated with 
galaxies. 

The most important strength of the HOD formulation of bias is its completeness. For a given 
cosmological model, the HOD tells us everything a theory of galaxy formation has to say about 
the statistics of galaxy clustering, in real space and redshift space, and on small, intermediate, 
and large scales [|. Conversely, if we can determine the HOD empirically, we will learn everything 
that observed galaxy clustering has to tell us about the physics of galaxy formation. Moreover, 
the HOD provides a physically informative basis for interpreting discrepancies between predicted 
and observed galaxy clustering, or between predictions of different galaxy formation theories 
themselves. It would be more illuminating to learn, for example, that a given theory predicts too 
many red galaxies in halos of mass 10^^ — IO^^Mq than to learn that it predicts the wrong 3-point 
correlation function of such galaxies. Finally, since the HOD describes bias at the level of systems 
near dynamical equilibrium, empirical determinations of the HOD can take advantage of mass 
estimation methods that are inapplicable on large scales. For example, traditional virial methods 
and X-ray mass estimates of clusters can provide fairly direct constraints on P{N\M) at high M 
(see § 5). 

The halo occupation framework has a long history, initially in analytic models that described 
galaxy clustering as a superposition of randomly distributed clusters with specified profiles and a 



We discuss some caveats to this assertion at the end of § 2 
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range of masses (INeyman fc Scott 1952| ; [McClelland fc Silk 1977| ; [Peebles 197^ ). The expl osion of 
recent activity in this field has been fueled partly by the recognition that a combination of this 
approach with recently developed tools for predicting the spatial clustering of halos ( |Mo fc White 



1996) provides a powerful formalism for analytic calculations of dark matter clustering, which can 



be naturally extended to biased galaxy populations (e.g., ^eljak 2000| ; Ma fc Fry 200C ; Peacock 
& Smith 2000; [Scoccimarro et al. 2000| ; White 200 1[ ; and numerous other papers referred to in 
subsequent sections). Our own interest was sparked largely by the paper of Benson et al.[ (|2000 



see Kauffmann, Nusser, & Steinmetz 1997 and Kauffmann et al. 1999 for similar analyses), who 
demonstrated that they could reproduce the clustering of galaxies in their semi-analytic models by 
populating N-body halos according to a predicted P{N\M). Furthermore, they showed that the 
predicted clustering depends not only on the complex mass dependence of the mean occupation, 
but also on finer details of sub-Poisson fiuctuations about the mean. This result illustrates the 
power of the HOD to test detailed predictions of galaxy formation theories. Models of P{N\M) 
based on semi-analytic calculations of Benson et al. (2000| ) and Kauffmann et al. (1999| ) have 
been incorporated into several of the papers cited above, and some recent papers have presented 
predictions of hydrodynamic simulations for P{N\M) ( [White, Hernquist, fc Springel 2001 ; 



Yoshikawa et al. 2001). We will compare predictions of P{N\M) from hydrodynamic simulations 



and semi-analytic calculations in a future study (Berlind et al., in preparation). 

The HOD description can be contrasted with another widely used approach that characterizes 
bias in terms of the correlation between galaxy density and properties of the large-scale 
environment, such as mass density, temperature, and geometry. This "environmental bias" 
approach has been used to study the effects of generic biasing models on galaxy clustering 
statistics ( Weinberg 1995 ; [Mann, Peacock, fc Heavens 1998t [Dekel Sz Lahav 1999[ ; Narayanan, 
Berlind, & Weinberg 2000; Berlind, Narayanan, & Weinberg 2001) and to encapsulate predictions 
of hydrodynamic simulations and semi- analytic galaxy formation models ( Blanton et al. 199S| ; 



Cen & Ostriker 2000; Somerville et al. 2001; Yoshikawa et al. 2001). It has also led to valuable 



analytic results concerning the shape of the galaxy power spectrum and the influence of bias on 



higher-order clustering on large scales (Poles 1993 ; Fry & Gaztaiiaga 1993; Fry 1994; Juszkiewicz 
et al. 19951 ; [Scherrer fc Weinberg 199^ ; [Coles, Melott, fc Munshi 1999|) . However, this formulation 



cannot effectively describe bias on scales smaller than the smoothing length used to define the 
environment, and the choice of smoothing scale is, itself, rather arbitrary. In the HOD framework, 
there is some range of reasonable methods for defining halos, but the choice of p/p ~ 200 for 
a typical halo boundary is well motivated by the division between the infall and dynamical 
equilibrium regimes. Also, as already noted, this choice allows use of virial mass estimates to 
constrain P{N\M) empirically, while the large-scale matter density, which plays a fundamental 
role in environmental formulations of bias, is generally unobservable. 

Our methodology in this paper is to define generic models of the HOD, apply them to an 
N-body simulation of the inflationary cold dark matter scenario (see § 2), and investigate the 
dependence of galaxy clustering statistics on the HOD parameters. This numerical approach 
complements earlier analytic work on halo bias by considering a wider range of clustering statistics 
and HOD models, some of them not readily amenable to analytic calculations. We will interpret 
our results in light of the analytic formalism developed in other papers and in terms of some 
heuristic analytic models presented here as a guide for understanding. The next section deflnes our 



- 4 - 



HOD prescriptions more formally. We then devote considerable attention, in § 3, to the two-point 
correlation function £,{r), because of its intrinsic importance and in order to illustrate some general 
features of the way the HOD influences galaxy clustering. We consider other clustering statistics 
in § 4, and properties of galaxy groups in § 5. A general theme emerging from these results is that 
different statistics respond to different features of the HOD, implying that precise measurements 
of galaxy clustering can in principle yield an empirical determination of the HOD. Achieving this 
goal in practice requires a scheme for getting a good first guess at the HOD that will reproduce 
the observed galaxy clustering. We outline such a scheme and present an illustrative test in § 6. In 
its present form, this scheme assumes that the underlying cosmology is known, but we speculate 
on prospects for breaking the degeneracy between bias and cosmology. We summarize our results 
in 5 7. 



2. Models of the Halo Occupation Distribution 

In the HOD framework, the relation between the galaxy and matter distributions is fully 
defined by 

(1) the probability distribution P{N\M) that a halo of virial mass M contains N galaxies, 

(2) the relation between the galaxy and dark matter spatial distributions within halos, and 

(3) the relation between the galaxy and dark matter velocity distributions within halos. 

We use the term "Halo Occupation Distribution" (or HOD) to refer to all three of these aspects. 
Each individual class of galaxies (defined, for example, by luminosity and color ranges or by 
morphological type) has its own HOD. 

For this study we have used the high resolution GIF N-body simulations carried out by the 
Virgo consortium (Jenkins et al. 199§|). The particular model we have focused on is the flat 



ACDM model with = 0.3, = 0.7, Hq = 70 km s Mpc , and a spectral shape parameter 



F = 0.21 (in the parameterization of Efstathiou, Bond, & White 1992). The simulation follows the 



evolution of 256^ particles, each of mass 1.4 x 10^^ h ^Mq, in a comoving box of size 141.3/i^^Mpc. 
The rms mass fluctuations on a scale of 8/i~^Mpc are erg = 0.9, in agreement with the observed 
abundance of clusters ( |Eke, Cole, fc Frenk 199^ ). Gravitational forces are significantly softer than 
1/r^ on scales ^ 30/i~^kpc. 

We identify halos in the dark matter distribution using a "friends-of-friends" (FoF) algorithm 



( Davis et al. 1985 ) with a linking length of 0.2 times the mean inter-particle separation, and we 
only consider halos consisting of 10 or more particles. This means that the smallest halos we can 
resolve have a mass of 1.4 x W^^h~^MQ. 

For our present purposes, an HOD prescription amounts to a recipe for selecting a subset 
of the dark matter particles in these halos to represent the galaxy population of the simulation. 
Most of the models we show in this paper are tuned to produce a galaxy population with a space 
density of fig = O.Ol/i^Mpc"'^. This space density corresponds to galaxies of luminosity L ^ 0.5L^, 
assuming no further cuts in color or morphology. Our HOD prescriptions consist of the following 
features: 

1- NavgiM) - We consider two types of models for this relation, power laws and broken power 
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laws. In the power-law models, the mean number of galaxies that populate dark matter halos of 
mass M is 

^ _ / if M < M^^^ 

^^avg j (M/Mi)" otherwise, ^ > 

where a is the power-law index, Mmin is the cutoff halo mass below which halos cannot contain 
galaxies, and Mi sets the amplitude of the relation and corresponds to the mass of halos that 
contain, on average, one galaxy. In the broken power-law models, 

if M < Mmin 

iVavg = { {M/MiY if AVin <M< Merit (2) 

{M/M[Y otherwise, 

where a and (3 are the low and high mass power-law indices, M^rit is the halo mass at which the 
power-law slope breaks, and M[ is required by continuity to be M[ = M'f^^'^M^-^^°'^^\ For a 
given HOD model, the value of Mi is chosen to produce a galaxy population of the desired space 
density. 

The parameters Mmin, «! /3j and M^it are directly related to the efficiency of galaxy formation 
as a function of halo mass. For example, bright galaxies cannot form in halos below a certain mass 
because these halos do not contain enough cold gas, hence the need for Mmin- The simplest form 
of iVavg(M) would have N proportional to M (a = 1) for M > Mmin. However, many physical 
mechanisms can alter the efficiency of galaxy formation as a function of halo mass. For example, 
the typical cooling time for gas increases with halo mass, and this suggests that a < 1. This effect 
could be counteracted by an earlier average formation time of galaxies that end up in high mass 
halos. Galaxy mergers may alter galaxy numbers preferentially in intermediate mass halos, where 
mergers are most frequent. Any given physical process can affect A''avg(M) differently for different 
galaxy classes. For example, mergers decrease the number of low luminosity galaxies but increase 
the number of high luminosity galaxies. Morphological transformations, likewise, increase the 
numbers of one galaxy type while decreasing those of another. We will, therefore, be able to learn 
much about the processes that determine galaxy properties by comparing the HOD for different 
galaxy classes. 

2. P(A^|A''avg) - Once A'^avg is determined, the actual number of galaxies that occupies any 
given halo is drawn from a probability distribution P(A^|A''avg). We consider three such probability 
distributions: (a) a Poisson distribution; (b) a very narrow distribution, which we call "Average", 
where the actual number of galaxies is the integer either above or below N^^-yg, with relative 
frequencies needed to give the required mean (this is identical to the "Average" distribution used 
by Benson et al. 200C| ); (c) a Negative Binomial distribution, which is substantially wider than 



Poisson. Although it is possible for the form of P(A'^|A''avg) to vary with M, we do not consider 
such models in this paper. Figure Q shows P(A'^|M) for a particular HOD prescription (power-law 
^sLvg{M) with Minin = 2.8 X IO^^H^^Mq, a = 0.5). Each point represents the number of galaxies 
chosen to occupy a particular halo in the dark matter distribution. The two panels show the 
difference between assuming an Average and a Poisson i-'(A^|A'avg)- 

In the high halo mass regime, P(A^|A''avg) should depend on the statistics of halo merger 
histories ( [Lacey &: Cole 1993| ). The distribution of merger histories for halos of a given mass 
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Fig. 1. — P{N\M) relation for two particular HOD models. Each point represents the number 
of galaxies that occupy a single halo in the N-body simulation. Points for halos that contain no 

galaxies arc arbitrarily placed at logA^ = —0.5. The HOD prescriptions shown have a power-law 
iVavg(M) with a = 0.5 and Mmin = 2.8 x lO^/i-^M© (see definitions in § 2), as indicated by the 
solid lines. The P(iV|iVavg) distribution are Average (top panel) and Poisson (bottom panel). 
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will produce a resulting distribution of galaxy numbers for those halos. In the low mass regime, 
each halo is expected to contain only one galaxy, but whether that galaxy would pass a given 
luminosity threshold depends on the gas cooling and star formation history of that halo. Therefore, 
P(A'^lA^avg) depends on the regularity of galaxy formation in halos of a given final mass. 

3. Central galaxy - Once the actual number of galaxies N that occupy each halo is determined, 
we must specify how these galaxies are distributed within halos. The first step in this process is 
to specify whether or not there must be a galaxy at the center of each halo for which N > 0. We 
investigate both these cases. If we force a galaxy to sit at the halo center, we place it at the halo 
center of mass and assign it the mean halo velocity. 

4. Galaxy concentration - We allow for the possibility that galaxies are more or less spatially 
concentrated than the dark matter within halos. We implement such models by selecting "galaxy" 
particles with probability P oc r'^^, so that, on average, 

p,{r)/pUr)=r''^. (3) 

This prescription imposes a difference A7 in the logarithmic slopes of galaxy and dark matter 
profiles without imposing any specific form or symmetry on the galaxy distribution in the halos; 
the galaxy distribution will inherit the geometry of the dark matter. A non-zero A7 can be 
applied together with the central galaxy prescription or on its own. 

5. Velocity bias - Finally, we allow for the possibility of velocity bias within halos. The mean 
velocity of galaxies in a halo should not differ from that of the dark matter, since both components 
are responding to the same large-scale gravitational field. However, the galaxies in a halo might 
have a higher or lower velocity dispersion than the dark matter particles at the same locations. 
We define a velocity bias factor ay through the relation 

Vg = Vh + a^(vm - Vh), (4) 

where Vg and Vm are the velocities of the galaxies and the dark matter particles that they are 
assigned to and is the mean center-of-mass halo velocity. For example, if a^, = all the galaxies 
within a halo have the mean halo velocity, while if a^, = 1 galaxy velocities trace the dark matter 
velocities. Our central galaxy and galaxy concentration prescriptions impose some degree of 
velocity bias even if a^, = 1, the first because the central galaxy is assumed to move at v^, and the 
second because the typical dark matter velocities depend on radius if the halo is not isothermal. 

Prescriptions 3, 4, and 5 allow us to represent a number of physical processes that could 
affect the galaxy distribution in important ways. For example, dynamical friction could cause 
galaxies to sink to the center of a halo and end up with a spatial distribution that is more centrally 
concentrated and a velocity distribution that is colder than that of dark matter. In addition, if 
galaxies form near the centers of their original parent halos, they can inherit spatial and velocity 
bias as a result of incomplete relaxation during the merging of these halos into a larger common 
halo ( [Evrard 1987 ). On the other hand, galaxy mergers happening at the centers of massive 



halos could reduce galaxy numbers in those regions, thus causing galaxies to be less centrally 
concentrated than dark matter (A7 > 0). Also, an a^, > 1 velocity bias could arise as a result of 
preferential destruction or merging of galaxies that have a low velocity. 
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Following the above steps, we have created a large number of galaxy distributions spanning a 
wide space of HOD parameters. All of these galaxy distributions come from the same dark halo 
population and differ only in the HOD. We have calculated a variety of clustering statistics for 
each of these galaxy distributions in order to test the sensitivity of each statistic to features of 
the HOD. It would be impractical to show all statistics for all of our HOD models, so in each of 
the following sections we focus on a subset of models that illustrate the sensitivity of the statistic 
under examination. Our approach of populating N-body simulations according to an HOD with 
power-law A^avg(-^) is similar to that used by |Jing, Mo, fc Borner (1998 ), Jing fc Borner (199^ ), 



and Jing, Borner &: Suto (200^ ) in their modeling of measurements from the LCRS and PSCz 



redshift surveys. However, the implementations are different, and our HOD model is substantially 
more general: we allow variations of Mmin and the form of P(A^|A'^avg) in addition to variations in 
a, and we allow the possibility of spatial and velocity biases within halos. 

We asserted in § 1 that the HOD formulation can provide a complete statistical description 
of the bias between galaxies and mass. Underlying this assertion is an implicit assumption that 
the galaxy content of a halo of virial mass M is statistically independent of that halo's larger 
scale environment. This assumption is supported by the N-body simulation results of [Lemson fc 



Kauffmann (1999), who find, in agreement with predictions of the excursion set model ( [Bond et 



al. 1991 ), that halos of fixed mass in different environments have similar properties and formation 
histories, although the halo mass function is itself systematically shifted towards higher mass 
halos in high density regions. However, the alleged independence of halo histories and large 
scale environment merits more detailed theoretical investigation. The hypothesis that the HOD 
formulation is complete will ultimately be tested empirically, by seeing whether an HOD model 
can reproduce all facets of observed galaxy clustering when applied to a cosmological model 
consistent with other observational data. Additional implicit assumptions, that all galaxies reside 
in virialized dark matter halos and that the halo population itself is minimally affected by baryonic 
physics, seem well justified, though the latter deserves more thorough testing with hydrodynamic 
simulations 0. The specific parameterizations adopted here may not capture all of the important 
features of the true HOD, though they are flexible enough to produce a wide range of results and 
can easily be generalized to include, e.g., mass dependence of P(A^|A''avg) or A7. 



3. The Galaxy Correlation Function 



We begin our analysis with the two-point correlation function (,g{r), which plays a fundamental 
role in understanding galaxy clustering because it has been thoroughly studied as a function 
of galaxy type, color, and luminosity (see Norberg et al. 2001, 2002, Zehavi et al. 2001, and 



numerous references therein) and because its observed form is remarkably simple. For typical 
optically selected samples, ^g(r) is a power-law (r/rg)"'^ for separations r ^ 5/i~^Mpc, with 
ro ~ 5 — 6/i~^Mpc and 7 ~ 1.8. More luminous galaxies have a larger rg and similar 7, while 



^The second assumption would not hold if we defined halos at a much higher overdensity, since dissipative 
collapse of baryon clumps increases their ability to retain surrounding dark matter concentrations within groups 
and clusters. However, from our point of view, high density halos surrounding individual galaxies are substructure 
within overdensity 200 halos, so they are described statistically by the HOD itself. 
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redder or early-type galaxies have a larger 7 and a higher clustering amplitude on small scales. 



Cosmological N-body simulations show that CDM models do not predict a power-law matter 
correlation function (e.g., Jenkins et al. 199^ ), and analytic theory ( Hamilton et al. 1991| ; Peacock 
fc Dodds 1996| ) implies that the linear theory power spectrum would have to contain a specially 
tuned feature in order to yield a power-law ^(r) on non-linear scales. If the primordial power 
spectrum is a smooth function, as expected on theoretical grounds, it appears that scale-dependent 
bias must transform the curved matter correlation function into a power-law galaxy correlation 
function. Remarkably, galaxy distributions predicted by hydrodynamic simulations, by "sub-halo" 
analyses of high resolution N-body simulations, and by semi-analytic models applied to N-body 
halos all yield power-law galaxy correlation functions, at least for some reasonable choices of 
cosmology and galaxy definition parameters ( [Pearce et al. 199£ ; Dave et al. 2000 ; Cen fc Ostriker 
20001 ; [Yoshikawa et al. 2001| ; |Colm et al. 19991 ; [Kauffmann et al. 1999| ; [Benson et al. 200^ ; 
Somerville et al. 2001| ). 



In the context of HOD bias, we would like to know whether a power-law (,g{r) follows 
from some simple and generic property of the HOD, such as a particular galaxy profile in high 
multiplicity halos, or whether it demands finely tuned parameter choices. We also want to know 
more generally how the amplitude and shape of ^g(r) depend on parameters of the HOD. 



3.1. Analytic Discussion 



Recent papers ( geljak 200d| ; |Ma fc Fry 2000| ; gcoccimarro et al. 2000] ; gheth et al. 20011) 



present a fairly complete analytic theory of the galaxy correlation function in the halo bias model. 
The basic idea is to add the "1-halo" term representing galaxy pairs within a single halo to the 
"2-halo" term representing pairs in separate, spatially correlated halos (as done in a different 



context by pcherrer fc Bertschinger 199l| ). The full analytic theory becomes simple in the Fourier 
domain, where convolutions of the halo profile transform into multiplications. As a guide to 
interpreting our numerical results, we begin with a complementary, more approximate discussion 
of correlations in real space. 

For separations larger than the virial diameter of the largest halos, all pairs must come from 
galaxies in separate halos. Mo Sz White (199^ ) derived an analytic approximation for the bias 



factor of halos hf^{M) as a function of halo mass M using the Press-Schechter (1974| ) formalism 



Above the characteristic mass in the Press-Schechter mass function, halo formation is enhanced 
in regions of high background density (and suppressed in underdense regions), so hh{M) exceeds 
unity and increases rapidly with increasing mass. Halos with M < are weakly anti-biased 
because they merge into more massive systems in overdense regions. [Jing's (1998| ) numerical 
results yield bh{M) 0.7 — 0.8 for M <C M*. We assume that the cross-correlation between halos 
of mass Ml and M2 is 

^i2{r)=bh{Mi)bh{M2)Uir), (5) 

where ^mif) is the mass correlation function. 

Consider a halo of mass Mi . The mean number of excess galaxies in a volume dV at distance 
r from the halo is obtained by integrating over the differential halo mass function dn/dM (which 
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has units of number density per unit mass), weighting each halo by the mean number of galaxies 
Ng^-vgiM) and by the biased correlation factor: 

iVexcess = / dM^-— N^^g{M2)bh{Mi)bh{M2)U{r)dV. (6) 
-'M,„i„ dM2 

If the number density of galaxies is fig, then the number density of correlated galaxy pairs is 

1 1 dn 

-nlUr)dV = -/.. dMi-^iV^vg(Mi)iVe: 



2 7m^,„ dMi 

2^M„,i„ dMi Jm^,^ dM2 
1 /■°° ^ dn ^ . r°° ^ dn 



dMi-^N^^^{Mi) / dM2-—N^.^{M2)bh{M{)hh{M2)U{r)dV 



-U{r)dV / dMi— -iVa,g(Mi)6/,(Mi) / dMa— -iVavg(M2)6ft(M2)(7) 

^ JM ■ dM-i JM ■ (llVlo 

where the factor of 1/2 corrects for double counting of each pair. Thus, at large separations we 
expect 



igir) = b'Uir), b = n-' j dM—N,^^{M)bh{M), (8) 



dn 

'M^,^ dM 

and the galaxy bias is just the weighted value of the halo bias. 



On small scales, the correlation function is dominated by the 1-halo term representing galaxy 
pairs that reside in the same halo. The total number density of such pairs is 

dn(N(N-l)),. nl /-sfimax 
npair,ih = / rfM— ^ ^ '"^ ^ -f / ig{r)A7:r^dr, (9) 

JM,-nin dM I 2 Jo 

where {N{N - l))^^^ = dNP{N\M)N{N - 1) and the second equality would be exact if all 
correlated pairs out to the maximum halo virial diameter 2i?max came from the 1-halo contribution. 
( Bullock et al. 2002] apply a similar argument to Lyman-break galaxy clustering.) More generally, 
we can write 

2 Jo ^ ^ ' Jm^,^ dM 2 \2R^,J ' ^ ^ 

where ^ih('^) refers specifically to the correlation function associated with galaxy pairs in the same 
halo and the function F(r/2i?vir) represents the average fraction of galaxy pairs in halos of mass 
M that have separation less than r (a fraction r/2i?vir of the virial diameter). Since the only r 
dependence on the r.h.s. of equation ([l0|) appears in F(r/2i?vir), we can differentiate with respect 
to r to find 

' dM 2 2i?,ir(M) \2R^,J ^ ' 

While the large scale bias factor (eq. Q) depends only on A'^avg (Af ) , the correlation function 
in the 1-halo regime depends on the 2nd factorial moment {N{N — 1))^/, and hence on the 
form of P(A^|A''avg)- A Poisson distribution of mean {N) has variance (N^) - {Ny = {N), so 
(A^(A^ - 1))^ = {N'^)m - Wm = (^)m = ^avg(^)- Our Average distribution, on the other 
hand, has {N{N - l))j^^ = N^^^{MI) - iVavg(Af) (even when A^avg(M) is not an integer), which 
can be substantially lower than the Poisson value when A''avg(Af) is small. The small-scale galaxy 
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correlation function also depends on the halo profile through F(r/2i?vir) — more concentrated 
halos have a larger fraction of pairs at small r, boosting the small-scale correlations at the expense 
of slightly larger separations. 

Figure Qa shows the cumulative pair fraction F(r/2i?vir) for one of our HOD models. We take 
the virial radius of each halo to be 



/ 3M 



where pm is the mean mass density and M is the halo mass. We measure the separation 
distribution from the numerically realized galaxy distribution and average over all halos (solid 
line), as well as halos in narrower mass bins. With separations scaled to virial diameters, the 
function F(r/2i?vir) is nearly independent of halo mass, though there is a slight trend of higher 



concentration at lower mass, as expected from N-body studies of halo profiles (Navarro, Prenk, & 



White 1997; Bullock et al. 2001). All pairs have separations less than 2i?vir, and about half the 



pairs have separations less than i?vir/2. 

The halo-averaged form of F{r /2R^\^) should be virtually independent of P{N\M)^ since 
the cumulative pair distribution is insensitive to halo mass. However, F{r /2Rvh) is sensitive to 
the relative distribution of galaxies and dark matter within halos. Figure ^ shows -F(r/2i?vir) 
averaged over all halos for HOD models in which the galaxies are more centrally concentrated 
(A7 = —1; short-dashed curve) or less centrally concentrated (A7 = -|-1; long-dashed curve) than 
the dark matter. As expected, F{r /2R^\^) rises faster for more central concentration, since most 
pairs lie close to the halo center. The dotted curve in Figure |2|b shows a model with A7 = 
but a central galaxy in every halo (for which > 0), which yields a result intermediate between 
the A7 = and A7 = —1 models. In this case, unlike all the others, we expect -F(r/2i?vir) to 
depend strongly on halo mass, since the central galaxy contributes to a significant fraction of pairs 
in low multiplicity halos but not in high multiplicity halos. For central galaxy models, it is more 
informative to rewrite equation (|lO|) in the form 

^ r e,,(0 W^. = r dM^ \{{N-l){N-2)) /^X _ 

(13) 

where F(r/2i?vir) now refers only to pairs that do not involve the central galaxy and Fc{r /2R^i^), 
which is simply the cumulative halo profile itself, refers only to pairs that include the central 
galaxy. Here F{r /2R^\r) and Fc{r /2R^\r) should both be approximately independent of halo mass, 
but Fc rises faster than F, and the relative weight of the two terms depends on N . We see from 
equation ( [T^ ) that a central galaxy should be important only on scales where low multiplicity 
halos make a significant contribution to the correlation function. 



3.2. Numerical Results 

Figures |3[^ show correlation functions for galaxy distributions created using a variety of HOD 
models that all have a power-law N^^g{M) (as defined in eq. |]l|). Each figure demonstrates the 
dependence of ig{r) on a particular feature of the HOD by displaying models that differ only in 
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Fig. 2. — The cumulative fraction of galaxy pairs within halos as a function of separation, 
F(r/2i?vir)! as defined in equation ([To|). Panel (a) shows F(r/2i?vir) ^ov halos of different mass 
ranges in the case of an HOD model with a = 0.5, Mmin = 2.8 x lO^^h^^ Mq, and an Average 
P(A^|A'^a,vg)- Panel (b) shows F{r /2R^\r) averaged over the halo population for HOD models with 
the same P{N\M) as panel (a), but with variations of the galaxy spatial distribution within halos. 
There are three HOD prescriptions with different values of A7 and one where a galaxy is forced to 
lie at the center of each halo for which > 0. 
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that feature. We note that models with different values of Afmin or a also have different values of 
Ml, since this parameter is used to fix the mean number density of galaxies to fig = 0.01/i^Mpc~^. 
For purpose of comparison, each figure also shows the dark matter correlation function and the 
real space ^^(r) inferred by Baugh (1996| ) from angular clustering in the APM galaxy survey 



(iMaddox et al. 19900 - 



Figure ^ shows the effect on ^g(r) of varying M^ain- On large scales, only the amplitude of 
the correlation function is affected, with higher values of Mmin having a slightly larger amplitude. 
Since the number density of galaxies remains fixed, the result of increasing Mmin is to remove 
galaxies from low mass halos and redistribute them into halos of mass M > M^[^. In terms of 
equation (^), this means that weight is taken away from the lower mass halos that have a smaller 
halo bias factor b^, so the galaxy bias b increases. On scales r ^ l/i~^Mpc, both the shape and 
amplitude of the correlation function are very sensitive to Mmin, with higher values of Mmin 
producing a steeper slope and higher amplitude of ^g{r)- Figure Q shows the effect of varying a. 
Qualitatively, increasing a has an effect similar to that of increasing Mmin, since both changes 
boost the fraction of galaxies in high mass, positively biased halos. However, the changes to the 
shape of ^g{r) are different in detail. We will return to a discussion of these effects shortly. 

Figure |5| shows the effect of varying P(A''|A^avg) while keeping A^avg(-M) fixed. It is clear that 
the impact on large scales is negligible, while on small scales the amplitude of Cgi''") increases with 
increasing width of the P(A^|A''avg) distribution. This behavior is expected because the large-scale 



bias factor (eq. fg]) depends only on N^vgiM), while Cih{r) depends on {N{N — 1))^,^ (eq. [lll|), 
which is larger for wider distributions that have the same A^avg- The impact of P(A''|A^avg) grows 
at smaller scales, where low mass (and hence low multiplicity) halos contribute to 

Figure ^ shows the effect of varying the distribution of galaxies within halos, while keeping 
P{N\M) fixed. The plot shows four models: one where galaxies are more centrally concentrated 
than the dark matter within halos (A7 = —1), one where galaxies are less concentrated within 
halos (A7 = +1), one where galaxies trace the dark matter within halos (A7 = 0), and a model 
where a galaxy is forced to lie at the center of every halo for which > 0. The effects of varying 
the spatial distribution of galaxies within halos are confined to small scales. As expected, the 
model where galaxies are more centrally concentrated within halos has a correlation function that 
is amplified on very small scales and depressed on scales corresponding to the virial size of halos. 
The model where galaxies are less centrally concentrated exhibits the opposite behavior. Forcing 
halos to have central galaxies only affects ^g(r) by increasing it slightly on the smallest scales, 
as we expect from equation (p!3|). Relative to changes in P{N\M), the effect of changing the 
distribution of galaxies within halos is small, at least on the scales r > 0.1/i~^Mpc considered here. 



To better understand the behavior in Figures |^j6|, it is helpful to decompose S,ih{f) into 
contributions from halos in different mass ranges. Figure |^ presents such a decomposition for 
several HOD models (see figure 2 of peljak 2000 for a similar decomposition in the Fourier domain). 



Figure 0a shows a power-law Aavg(M) model with a = 0.5 and Mmin = 2.8 x 10^^ h~^MQ that 
has a Poisson P(A|Aavg). The thick solid curve shows the correlation function for this galaxy 
distribution (also seen in Figures ^-§), and the thin solid curve represents the contribution from 
galaxy pairs within the same halo (the 1-halo term). As expected, the 1-halo term dominates on 
scales up to the virial size of typical halos and drops off quickly at larger scales. The remaining 
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Fig. 3. — Influence of Mmin on the galaxy correlation function. Curves show galaxy correlation 
functions for HOD models with a power-law Na.-yg{M), a = 0.5, Average P(A^|A''avg), and different 
values of Mmin, which are listed in the legend in units of W^^h~^ Mq. The solid curve shows the 
mass correlation function, and points show the correlation function measured from the APM galaxy 
survey ( [Baugh 1996 ). 
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Fig. 4. — Influence of a on the galaxy correlation function. Curves show galaxy correlation functions 
for HOD models with a power-law iVavg(M), Mmi„ = 2.8 x W^^h'^MQ, Average P{N\Ns_yg), and 
different values of a, which are listed in the legend. 
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r [h~^Mpc] 

Fig. 5. — Influence of -P(iV| A'avg) on the galaxy correlation function. Curves show galaxy correlation 
functions for HOD models with a power- law Afavg(M), Mmin = 2.8 X lO^i/j-^Mo, a = 0.5, and 
different forms of P(A'^|A'^a,vg)) which are listed in the legend. 
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Fig. 6. — Influence of the galaxy profiles within halos on the galaxy correlation function. As in 
Figs. the dotted curve shows ^g{r) for a model with Mmin = 2.8 x W^^h~^MQ, a = 0.5, 
Average P(A^|A^avg)) and galaxies tracing dark matter within halos. Short-dashed and long-dashed 
curves show results for models in which galaxy profiles are respectively more concentrated or less 
concentrated than dark matter profiles (A7 = — 1 or A7 = +1, see eq. Q). The dot-dashed curve 
shows a model in which the first galaxy of each halo lies at the halo center and subsequent galaxies 
have the same profile as the dark matter. 
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four curves show the contribution to the 1-halo term from halos with M = 10 — 10 h Mq 
(dot-dashed curve), 10^^ - lO^^/i-^M© (dotted curve), 10^^ - lO^^/i-^M© (short dashed curve), 
and 10^^ — 10^^ /i~^Mq (long dashed curve). Each curve is highest at the smallest scales and drops 
off at larger scales. However, the curves for high mass halos start at lower amplitude and extend 
to larger r, since pairs are spread over a larger virial volume. Consequently, galaxy pairs in low 
mass halos dominate Cgi^) &t very small scales, while galaxy pairs in high mass halos dominate 
^g(r) on scales corresponding to their virial radii. 

Figure ^ shows the effect of doubling Mmin while keeping all other HOD parameters fixed. 
Since fig remains constant, galaxies that were previously in low mass halos are redistributed to 
halos above the new value of M^i^. In Figure]^, the contributions to ^g{r) from galaxy pairs in 
the three highest halo mass bins uniformly increase. The lowest mass bin (10^^ — 10^^/i^^Mq) 
behaves differently because M^[^ lies within it, causing some of the halos in this bin to lose pairs 
and others to gain pairs, resulting in no net change to the contribution from that bin. Doubling 
Minin slightly increases the large scale bias factor, but the impact on ^ih('") is much greater, so the 
overall effect is to steepen ^g(r). 

Figure |^ shows the effect of increasing a while keeping all other HOD parameters fixed. The 
contribution of the high mass halos increases dramatically, whereas the contribution of the lowest 
mass halos drops. We can thus understand why increasing a has a "smoother" effect on ^g{r) than 
increasing Mmiru producing less distortion in the overall shape. First, increasing a has a larger 
impact on the large scale bias factor, so the 1-halo and 2-halo terms rise more nearly in step. 
Second, increasing a flattens the shape of Cih('') as it increases its amplitude, by redistributing 
pairs to halos with larger virial radii. The shape of ^ih('^) therefore stays closer to the extrapolated 
shape of Cgi^) from large scales. 

Figure]^ shows the effect of changing the P(A^|A''avg) distribution from Poisson to Average. 
The narrower distribution yields a smaller value of {N{N — 1))^^,/ for the same Ni^-vg(M), especially 
in low multiplicity halos. As a result, the 1-halo contribution to £,g{r) drops dramatically at small 
scales, r ~ 0.1 - 0.4/i-iMpc. The contribution of halos with M = 10^^ - IO^^H-'^Mq disappears 
completely because Mi exceeds 10^'^h~^ Mq, so with an Average P(A^|A^avg) no halo in this mass 
range can have more than one galaxy. The suppression of pairs in low mass halos allows £,g(r) to 
continue as a power law down to r = 0.1/i^^Mpc. 

The last two panels of Figure ^ show the effect of changing the spatial distribution of galaxies 
within halos. The HOD models have the same P{N\M) as the model shown in Figure so the 
number of pairs in halos of each mass range is the same as before. However, the radial separations 
of these pairs are squeezed towards smaller r for A7 = — 1 (Fig. |7^) and stretched towards larger 
r for A7 = -|-1 (Fig. |^), making the correlation function steeper or shallower at small scales. The 
shapes of these curves are directly related to the function F{r /2R^\^), defined in equation (|lO|) and 
plotted in Figure ^, which depends on the radial profile of galaxies within halos. 
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Fig. 7. — Contributions to the galaxy correlation function from pairs within halos of different mass 
ranges. Each panel represents a particular HOD model. The models are specified at the top of each 
panel, where Mmin is given in units of 10^^h~^MQ, and Poi and Ave represent Poisson and Average 
P(A^|A''avg) distributions, respectively. Each panel shows the full galaxy correlation function (thick 
solid curve), the correlation function including only galaxy pairs that lie within the same halos 
(thin solid curve), and the contribution to ^g{r) from pairs that lie only within halos of a certain 
mass range (remaining 4 curves). The mass ranges are indicated in panel (d). Also shown, for 
comparison, is the APM galaxy correlation function (squares). 
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3.3. Understanding the Observed Correlation Function 

In light of these results, what do we make of the observed power-law form of ^g{r)l One 
obvious conclusion is that a power-law ^,g{r) is not a generic prediction of HOD models applied to 
a ACDM cosmology; most of the models in Figures ^[]6| show clear departures from a power law. 
We can understand this behavior by considering Figure |^ and the analytic discussion in § 3.1. 
On large scales, the shape of £,g{r) is the same as the shape of £,m{T), and the amplitude of S,g{r) 
is determined by the imi'i') amplitude and the bias factor b (eq. On small scales, typically 
r ^ 0.5/i~^Mpc (Fig. |^), the 1-halo term dominates. In this regime, the shape and amplitude of 
ig{r) are governed by equation (|Tl|), which involves {N{N — '^)) rather than Ng^^^{M) and does 
not involve ^m(^^) explicitly at all (though and dn/dM are connected indirectly). Achieving a 
power-law ig{r) requires that the amplitude of i\h{r) place it on the extension of b'^Cmir), and it 
requires that the distribution of pair counts as a function of halo mass yield a power law of the 
same slope in the 1-halo regime. 

Suppose we have a model that achieves this somewhat delicate balancing act, such as the 
a = 0.5, Mmin = 2.8 x lO^^h'^MQ, Average model. Changing P(A^|A''avg) has no effect at large 
scales, since b depends only on A'^avg(-^), but it changes the amplitude and shape of Cih{f), 
destroying the power-law behavior. Changing Mmim with a compensating change in Mi to keep 
fig fixed, has different effects in the 1-halo and 2-halo regimes. If Mmin <^ M*, then the impact 
on b will generally be modest. However, raising M^[^ redistributes pairs that were previously in 
low mass halos to high mass halos, substantially increasing the 1-halo contribution of these halos 
and destroying the previous balance. Changing a has a more complicated effect than changing 
-^min because it can significantly alter b and because it changes both the amplitude and shape of 
^ih(r). However, these three changes generally do not work in a way that preserves a power law. 
Finally, Figure |6| shows that changes to the halo profile have relatively little leverage for modifying 
a non-power-law £,g{r) into a power law. Even the rather drastic profile changes represented by 
A7 = ±1 only alter £,g{r) on fairly small scales. A central galaxy restriction affects still smaller 
scales, since central galaxy pairs only matter in low-A^ halos with small virial radii. Central 
galaxies probably have an important influence on (,g{r) at r < 100/i^^kpc, but they are not 
fundamental to understanding the 0.1 — 0.5/i^^Mpc regime. More generally, adopting the "right" 
halo profile can never compensate for having the "wrong" P(N\M). 

To roughly quantify the "difficulty" of obtaining a power-law S,g{r), we have carried out 
a systematic survey of HOD models with power-law Na_-yg{M), Average or Poisson P(A^|A'avg), 
and Ml chosen in all cases to yield fig = O.Ol/i^Mpc"^. Fi gure P summarizes the results 
for a representative sample of models. Using an unweighted least-squares fit in the range 
0.1/i^^Mpc < r < 10/i^^Mpc, we determine the best-fit power-law parameters rg and 7 for each 
model correlation function. As expected, the slope 7 steepens with increasing Mmin or a, since 
these changes amplify the 1-halo term relative to the 2-halo term. The infiuence of Mmin is 
stronger in Average models than in Poisson models (compare Figs. |8|e and ^f) because {N{N — 1)) 
cuts off more sharply near ik/j^jj^ in an Average model. The correlation length tq generally increases 
with increasing Mmin or a (Figs. ^ and ^). However, the dependence on Mmin is weak when 
Mmin is small, because the bias factor integral (eq. |^]) remains dominated by low mass halos. 
The dependence becomes stronger when Mmin becomes large enough to significantly increase the 
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fraction of galaxies in halos with M > (and thus > 1). Increasing a, by contrast, always 
increases the relative numbers of galaxies in more biased halos, and the dependence of rp on a is 
fairly steady. 

The top panels of Figure |8| show the of the power-law fit. Since we used an unweighted fit, 
we scale the values of by a constant factor so that models with visually acceptable power-law 
correlation functions have < 1- (If we were trying to match specific observational data with 
known error bars, we would follow a more rigorous methodology, but that is not our purpose 
here.) As one would guess from the sampling of models shown in Figures only a small region 
of this HOD parameter space — Average models with Mmin ;^ 3 x 10^^/i~^Mpc and a ~ 0.4 — 0.6 
— yields power-law correlation functions. In more general terms, achieving a power-law ig{r) 
requires that a large fraction of galaxies reside in low mass (and low multiplicity) halos; otherwise 
the 1-halo portion of ig{r) is too high relative to the 2-halo portion. Average models are more 
successful than Poisson models because they suppress {N{N — 1)) in low mass halos and thus 
prevent a steepening of Cgi^) small r. Values of Mmin smaller than our resolution threshold 
(1.4 X 10^^ Mq) might allow a Poisson model to work, but we are modeling the population of 
galaxies with L > 0.5L*, and a lower value of Mmin would force some of these galaxies into very 
low mass halos, contradicting constraints from the Tully-Fisher (1977[ ) relation. In addition, the 



bottom panels of Figure |^ show that a small Afmin implies a large Mi (Mi > 3 x 10^^ /i ^Mq 
for Mmin < 2.8 X IO^^/j-^Mq). Since halos of mass ~ 2Mi still have a significant probability of 
hosting a single galaxy, a high Mi may force some low luminosity galaxies into fairly massive 
halos, again running afoul of Tully-Fisher constraints. 

This discussion suggests that the power-law Na.-yg{M) parameter space may itself be too 
restrictive, since it ties the fraction of galaxies in low mass halos directly to the relative distribution 
of galaxies among high mass halos (which affects the large scale bias and the larger separation end 
of the 1-halo regime). We have also examined the family of broken power-law models defined by 
equation (^), which (for a < (3) allow more galaxies in low mass halos for a given slope in the high 
mass regime. With a break point at M^it = 10^'^/i~^Mpc, we find acceptable correlation functions 
for some Average models with Mmin ~ 3 — 5 x IO^^H^^Mq, a ~ 0.2 — 0.3, f3 ~ 0.6 — 0.8, and for 
some Poisson models with Mmin ~ 1 - 2 x lO^^/i^^M©, a ~ 0.4 - 0.5, /? ~ 0.8 - 0.9. The solid 
line in Figure ^ shows the correlation function of one of these models, which matches the APM 
results quite well. Even with broken power-law iVavg(M), acceptable Poisson models require low 
values of Mmin- To the extent that such low values are implausible, we concur with arguments 
in previous papers ( [Benson et al. 2000| ; Seljak 20"00| ; peacock &: Smith 2000 ; Scoccimarro et al. 



2000) that sub- Poisson number fiuctuations in low multiplicity halos are important in explaining 
the observed form of ig{r). 

At scales r ^ 5/i~^Mpc, peculiar velocities severely distort the shape of the correlation 
function in redshift space. Observational evidence for a power-law ig{r) therefore rests on 
measurements of the angular correlation function w{9) or the projected redshift-space correlation 



function Wp{rp) (see, e.g., Baugh 1996| ; [Norberg et al. 2001 ; Zehavi et al. 2001 ). While a power-law 



^g(r) projects into a power-law w^O) or Wp{rp) ( [Limber 1954 ; Davis & Peebles 1983), one might 



worry that projection could wash out departures from a power-law S,g{r). We have checked that 
this is not the case by computing Wp{rp) for many of our HOD models. We find that any models 
that predict visually evident departures from a power-law (,g{r) also predict visually evident 
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Fig. 8. — Power-law fits to a sample of HOD model galaxy correlation functions. Each point 
corresponds to a specific P{N\M), where the x-axis shows M^am and the type of point represents 
a as shown in panel (b). The left-hand-side and right-hand-side panels contain models with an 
Average and Poisson P{N\Ng_yg), respectively. Panels (a) and (b) show values of from these 
fits. Since the fits are unweighted, the values of have been scaled so that models with visually 
acceptable power-law correlation functions have < 1- Panels (c) and (d) show the best-fit values 
of the correlation length vq. Panels (e) and (f) show the best-fit values of the slope 7. Panels (g) 
and (h) show the value of Mi for each model, chosen to yield a galaxy distribution of fixed number 
density fig = O.Ol/i^Mpc"^. 
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Fig. 9. — Real-space correlation functions of galaxy distributions for different HOD prescriptions. 

As indicated in the legend, there is one HOD model with a broken power-law NavgiM), and four 
power-law Na;yg{M) models with two values of the galaxy space density fig (see text for further 
discussion). Values of Mmin and Merit are given in units of 10^^h~^MQ, and values of fig are in 
units of h^Mpc~^. 
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departures from a power-law Wp{r.p). 

At higher luminosities, the observed ig{r) has a higher amplitude and a similar power-law 
slope. The higher amplitude is not surprising, since Mmin should be higher for more luminous 
galaxies. Retaining the power-law form of £,g{r), however, requires that other HOD parameters 
adjust in the right way. A higher luminosity threshold also implies a lower galaxy space density 
fig, so in general we expect both Mmin and Mi to be higher for more luminous samples. The 
influence of Mi on £,g{r) depends crucially on the form of the probability distribution P(A^[A'^a,vg)- 
For a Poisson model, raising Mi is equivalent to randomly sampling the galaxy population, which 
does not alter the correlation function. For an Average model, on the other hand, raising Mi 
reduces {N{N — 1)) by a factor larger than fig, thus reducing the 1-halo contributions to Cg{r). We 
demonstrate this point with an explicit example in Figure ^ The long-dashed and dot-dashed lines 
show £,g{r) for two models, one Average and one Poisson, with space densities fig = 0.01/i^Mpc~^. 
The two correlation functions show similar departures from a power-law at r 2/i~^Mpc. Raising 
Ml to reduce the number density to fig = O.OOl/i^Mpc"^ does not change the correlation function 
of the Poisson model (dotted line). However, the correlation function of the Average model drops 
dramatically on small scales and stays the same on large scales (short-dashed-line) , yielding a 
(noisy) power-law Cg(^)- To raise the amplitude of S,g{r) on large scales, one would want to increase 
-^min and/or a as well as Mi, changes that tend to boost the 1-halo term relative to the 2-halo 
term. The generic difficulty of finding HOD models that yield a high bias factor and a power-law 
^g(r) (see Fig. ^) suggests that suppression of 1-halo clustering by sub- Poisson fluctuations in 
low multiplicity halos is especially important for explaining the correlation function of luminous 
galaxies. 

The quantitative results in this Section (and the rest of the paper) apply to the specific ACDM 
cosmological model used for the N-body simulation described in § 2. The galaxy correlation 
function depends on the HOD and on the mass and spatial distributions of dark matter halos 
themselves. In terms of the analytic discussion in § 3.1, the cosmological model influences £,g{r) 
by determining dn/dM, bh{M), and ^m{r) (see eqs. ^ and |ll]). However, we expect that our 
conclusions about the influence of HOD parameters on £,g{r) and the generic requirements for 
obtaining a power-law correlation function would continue to hold for a fairly wide range of 
cosmological models. 

In HOD models, a power-law correlation function emerges from a balance of several competing 
effects and therefore requires rather specific combinations of parameters. This complex explanation 
of a simple result may at first seem unreasonably contrived, but we see no physically realistic 
alternative. A model with Ni^-vg{M) oc M above a mass threshold Mmin predicts a strongly 
scale-dependent bias in the non-linear regime if Mmin is large enough to exclude a significant 
fraction of the mass. Thus, even if the non-linear matter correlation function were a pure power 
law, that simplicity would be lost as soon as one demanded that luminous galaxies inhabit halos 
above some minimum mass. Eliminating the mass in low mass halos will in general change the 
1-halo and 2-halo regimes of the correlation function by different factors, so there is no reason 
to think that a power-law dark matter correlation function would produce a power-law galaxy 
correlation function if a significant fraction of the mass resided in halos that were too small 
to host galaxies above the luminosity threshold, even if the galaxies traced the mass in larger 
halos. Furthermore, obtaining a power-law £,rnif) requires features in the primordial matter power 
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spectrum that have no obvious physical motivation ( peacock fc Dodds 1996| ). 

We conclude instead that the observed power-law form of $.g{r) provides a strong constraint 
on the physics of galaxy formation, and that the success of semi-analytic models ( [Kauffmann et 



al. 1999; Benson et al. 2000), hydrodynamic simulations ( [Pearce et al. 199q ; |Cen fc Ostriker 2000|; 



Dave et al. 2000 ; Yoshikawa et al. 2001 ), and high resolution N-body simulations ( |Colm et al 



1999) in reproducing this form is entirely non-trivial. In a subsequent paper (Berlind et al., in 
preparation), we will show that semi-analytic calculations and SPH simulations of the ACDM 
model predict P{N\M) distributions that are similar to those of the broken power-law Average 
model shown in Figure which does indeed yield a power-law ^g{r). Improving observations of 
galaxy clustering already show evidence for statistically significant departures from a power-law 
correlation function on small scales at low ( Connolly et al. 200l| ) and high (porciani fc Giavalisco 
2001 ) redshift, in addition to the "shoulder" in S,g{r) at r ~ tq ( paztafiaga fc Juszkiewicz 2001 



and references therein). Such departures are theoretically expected at some level, and they may 
provide further important clues to the properties of the HOD. 



Other Statistics 



4.1. Galaxy-Mass Cross- Correlation Function 

Advances in wide-field digital imaging have recently begun to provide a new and direct probe 
of the relation between galaxies and dark matter, the galaxy-mass cross-correlation function 
^gm.{f)- This statistic can be measured using galaxy-gal axy lensing ([Brainerd, Blandford, &: Small 
19961 ; iGriffiths et al. 1996| ; |deir Antonio fc Tyson 199^ ; [Hudson et al. 1998| ; [Fischer et al. 2000| ; 



Smith et al.~200l| ) or by cross-correlating cosmic shear maps with galaxy light maps ( [Wilson, | 



Kaiser, & Luppino 2001). Since the lensing signal depends on the surface density in units of 
the critical density for lensing, rather than the cosmic mean density, these techniques measure a 
combination of Vtm and ^gm{f) rather than ^gmif) alone. Recent analyses of galaxy-galaxy lensing 
in the SDSS ( McKay et al. 2001 ) constrain ^gm{i^) out to r ~ l/i~^Mpc and demonstrate its 
dependence on galaxy luminosity, color, and environment. 

Theoretical calculations of £,gm{f) have been presented in the context of halo occupation 
models ([Seljak 2000| ), semi-analytic galaxy formation models ( puzik, &: Seljak 200 ip , and 
hydrodynamic simulations ( White, Hernquist, &: Springel 200l[ ). Much of our analytic discussion 
of ^,g{r) in § 3.1 carries over to ^gmii^), except that we are now considering the correlation between 
the galaxies and a second "population" of objects that trace the mass. On large scales, we expect 

Cgm{r) = hU{r), (14) 
with the same bias factor given in equation (|8|). On small scales, we expect a "1-halo" term 



analogous to equation (10), but with halos weighted by MN.^^g{M) instead of {N{N — l))jy/. An 



important consequence of this change is that the shape of P(A^|A''avg) should not influence igm{f) 
even in the 1-halo regime. 



Figure 10 shows the galaxy-mass cross-correlation function for a variety of HOD models, with 



each panel concentrating on a specific feature of the HOD. Figure 10a shows the effect of varying 
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-^min- On large scales, models with higher Mmin have a slightly higher bias factor, in agreement 
with our results for ^g(r) in Figure ^ On small scales, higher Mmin models have a significantly 
larger ^gmii^), since a larger fraction of their galaxies live in high mass halos. Increasing a raises 



igm{f) for the same reasons, as shown in Figure 10b. The influence of a extends over a wider range 



of scales, as it did for ^g(r) (see the discussion of Fig. |^ in § 3.2). Figure |lO| c confirms that igm{f) 
does not depend on the shape of P{N\Ns;^^), since the 1-halo and 2-halo contributions are both 
controlled by Ni^-^g{M). For A7 = — 1, galaxies are more concentrated in the central, high density 
regions of halos, boosting ^g-m{i~) on small scales (Fig. p!o|d). A central galaxy restriction has a 
similar effect, but the impact is confined to smaller r. Conversely, setting A7 = +1 suppresses 
^gm{f) by pushing galaxies to lower density regions of their parent halos. 

For r ^ 200/1- ^kpc, our approach of selecting N-body particles to represent galaxies 
necessarily underestimates igm{f) for a given HOD because it cannot account for the ability of 
condensed baryons to concentrate dark matter around them. One could do somewhat better by 
locating galaxies at density peaks of halo substructure, but this method still would not capture 
the effects of baryonic dissipation on the dark matter distribution. We expect our approach to 
be accurate at larger scales (note, for example, that the effect of a central galaxy restriction 
is negligible beyond 200/1" ^kpc) and to illustrate the correct qualitative dependence on HOD 
parameters even at small r. 



With observations of S,g{r) and (,gm{r), one can define a bias function b{r) = (,g{^)/t 



gm 



r . 



At least on large scales, we expect this ratio to equal the bias between the galaxy and matter 
correlation functions [^g{r) / (,m{'>^)]^^'^ , allowing ^mii^) to be inferred indirectly. Figure ^ confirms 
this expectation for several HOD models. In each panel, thick and thin solid curves represent the 
two bias functions for an HOD model with Mmin = 2.8 x 10^^h~^MQ, a = 0.5, and an Average 
P(A^|A''avg). Thick and thin dashed curves show corresponding results for a model with higher 
Mmin (O^-)' higher a (|lllb), or Poisson P(A^|A^avg) (^)- Ah of the curves go flat at r ;^ 4/i-iMpc, 
indicating that Cg(^) and S,gm{r) have the same shape as (,m{r). In this regime, the two definitions 
of b are in essentially perfect agreement. At smaller r, the bias functions are scale-dependent and 
do not track each other perfectly, though for a given model, they have similar shapes even in this 
regime. From Figure ^ we see that a similarity of shape between ^g{r) and Cgm{r) should not 
be taken as evidence that galaxies trace mass (Wilson et al. 2001). However, if measurements of 
igm{f) can be pushed to r ;^ 4/i-^Mpc, they should yield an accurate estimate of igir) / im.{f) , 
and even on smaller scales they can indicate whether galaxies are more or less clustered than the 
matter. Since the weak lensing signal depends on as well as £,gm{f)-, these measurements will 
generally constrain a combination of Jim and b rather than b alone. 
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Fig. 10. — Galaxy-mass cross-correlation functions for different HOD prescriptions. In each panel, 
the dotted curve shows Cgm{r) for a model with power- law iVavg(M), Mmin = 2.8 X IO^^H-^Mq, 
a = 0.5, Average P{N\Ns,vg), and A7 = 0. Other curves show Cgm{f) for models with different 
Mmin (panel a), different a (panel b), different P{N\Na,yg) (panel c), or different galaxy distributions 
within halos (panel d). 
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Fig. 11. — Effective bias functions defined by ratios of correlation functions for different HOD 
prescriptions. Thick curves represent {$,g/^m)^^^, where and ^rn are the galaxy and mass 
autocorrelation functions, respectively. Thin curves represent (Cg/^gm)^^^) where ^g„i is the galaxy- 
mass cross-correlation function. In each panel, solid curves represent an HOD model with power-law 
iVavg(M), Mmin = 2.8 X lO^^/i-^M©, a = 0.5, and Average P(A^|A^avg)- Dashed curves represent 
models with different Mmin (panel a), different a (panel b), or different P{N\Ng;vg) (panel c). 
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If one's goal is to determine the HOD given a known cosmology, then the principal virtue 
of ^gm{f) is its sensitivity to A''avg(M) in a regime where S,g{r) is sensitive to {N{N — 1))^/. By 
combining the two statistics, one can discriminate between models that have different P{N\Ns;yg), 
even if A^avg(-^) is adjusted so that the models produce similar galaxy correlation functions.^ 
However, independent data will not determine cosmological parameters perfectly, and the more 
crucial role of galaxy-galaxy lensing measurements will be to help break the degeneracy between 
the HOD and the background cosmology, by providing direct estimates of the masses of the halos 
that galaxies inhabit. We return to this point in § 6. 



4.2. Bispectrum 

Higher order clustering statistics complement the information in the two-point correlation 
functions ^g{r) and igm{f)- Because rare, high multiplicity halos make a larger contribution to 
high order statistics, we expect them to be more sensitive to the occupation of massive halos 
and to the high-iV tails of P(A^|A^avg)- |Ma &: Fry (200q ) and [Scoccimarro et al. (2000D give 



fairly comprehensive analytic discussions of the three-point correlation function and its Fourier 
transform, the bispectrum, in the halo occupation framework. The basic principles are similar to 
those that govern ^g{r) and Pg{k) — for example, small scale correlations are dominated by a 
1-halo term, in which halos of mass M contribute in proportion to the mean number of galaxy 
triples per halo, {N{N — 1)(A^ — 2))^^. However, the mathematics of three-point correlations is 
considerably more intricate, in part because there are three terms 2-, and 3-halo) contributing 
to the galaxy correlations, and in part because the bias of the halos is more complicated. On large 
scales, where second-order perturbation theory is valid, the relation between the halo and mass 
three-point correlations depends on two bias factors, one describing the ratio of rms halo and mass 
fluctuations and one describing the non-linearity of the relation between halo number and mass 
density contrast. 

To facilitate comparison with other work, we focus here on the reduced bispectrum, 

n(\. V 1. ^ - g(ki,k2,k3) I. _ IT. I 

g(ki, k2, k3) - + P{k,)P{k,) + P{k,)P{k,y - ^'^^ 

where i?(ki,k2,k3) is the bispectrum itself, and P{h) is the power spectrum. The three 
wave-vectors form a triangle ki -|- k2 + ks = 0, so they can also be specified by the magnitude 
ki, the ratio ki/k2, and the angle 9 between ki and k2. In second-order perturbation theory, the 
reduced matter bispectrum Qm is independent of the scale ki, but it depends on the triangle shape, 
with higher Qm at the elongated configurations {9 ~ 0° and 9 ~ 180°) favored by anisotropic 



gravitational collapse (Fry 1984). In bias models where the galaxy density contrast is a local 
function of the mass density contrast, 5g = f{6m), a second-order calculation yields 

= ^ + (16) 



^The constraint of fixed number density greatly restricts one's ability to trade ofi^ A'^avg(M) against P(A''|A'avg), 
since J dM ■j^Na.vgiM) must remain constant. However, fgm(r) provides more discriminating power because different 
scales constrain the behavior of A'avg(A'f) in difi^erent halo mass ranges. 
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where bi = /'(O) and 62 = ^/"(O) ( [Fry 1994 ; see also Fry &: Gaztanaga 1993 ; Juszkiewicz et 
al. 1995| ). On smaller scales, the value of Qm rises and then plateaus at a higher level, and the 
dependence of Qm on triangle shape decreases ( [Ma &: Fry 200d| ; Scoccimarro et al. 200q ). 



To compute the reduced bispectrum, we form continuous density fields by cloud-in-cell 
(CIC) binning the discrete mass or galaxy distributions onto a 200^ grid and calculate P{k) 
and i?(ki,k2,k3) using a Fast Fourier Transform (FFT). Figure |l2| shows Qeq{k), the reduced 
bispectrum as a function of scale in the equilateral triangle case {ki = ^2 = ^3)1 for a variety of HOD 



models and for the unbiased mass distribution. For the same set of HOD models, Figure 13 shows 
Q{9), the reduced bispectrum as a function of triangle shape, in the case ki = 2/c2 = l/iMpc~^. 
The fundamental mode of the 141.3/i~^Mpc simulation cube is kf = 0.0445/iMpc~^. In Figure 12 
Qcq{k) curves for the mass distribution begin at ki = 3kf, and in Figure 13 /ci = 2A;2 = 22.5kf. 



Figures 12a and 13a show models with a = 0.5, Average P(A^|A''avg), and varying values of 
Mjnin. Because of the limited size of the simulation volume, the mass Qeqik) does not show the 
low-A; plateau predicted by perturbation theory very clearly, though it does show the usual rise 



into the non-linear regime and plateau at high-/c. The scale ki adopted in Figure 13 is sufficiently 
non- linear that most of the shape dependence in Qm{(^) has disappeared, though there is still 
a slight enhancement for elongated triangles. The galaxies have much lower values of Qeq{k). 
In the context of equation (|l6[), we interpret the depression of Qeqik) as a sign of negative b2, 
since the galaxies are anti-hiased in an rms sense (compare the galaxy and mass correlation 
functions in Figure ^), and this anti-bias on its own would tend to boost Q (bi < 1). Turning 
to non-equilateral triangles in Figure 13 a, we see that Qg is strongly depressed relative to Qm for 
30° ^ G ^ 140° but rises towards Qm near 9 = 0° and, especially, near 6 = 180°. As a result, 
the shape dependence of Qg is stronger than that of Qm- We are unsure how to explain the 
shape dependence of the suppression of Qg, though a greater impact of bias on roughly equilateral 
configurations than on elongated configurations is not surprising. 

In Figure ^a, we see that changing Mmin has no significant impact on Qcq{k) in these 
models. While Beq{k) and P{k) both increase for larger Mmin, the two effects cancel in the ratio 
that defines Q. Figure p^ a shows a modest dependence of Q{9) on Mmin for 6 ^ 140°, perhaps 
because the signature of elongated structures is stronger when galaxies are distributed among 
lower mass halos (smaller Mmin). Figure [l2[ b shows the dependence of Qeq{k) on a; while Beq{k) 
and P{k) both increase for larger a, Beq{k) grows faster, and Qeqik) rises closer to the value of 
Qm- Because low a suppresses Q{9) primarily at intermediate values of 9, the triangle-shape 
dependence weakens as a increases (Fig. |l3|b). 



Figures 12c and 12d show the dependence of Qeq{k) on the form of P(A^|A''avg) for 
Mrain = 2.8 X IO^^H^^Mq and a = 0.5 and 0.8, respectively. The dependence is minimal in the 
first case but significant in the second, with the Average distribution producing the highest Qeqik) 
at high k. While the narrower width of the Average distribution leads to a smaller Beqik), it 



causes a still larger reduction in the denominator of equation (15), thus boosting Q. The greater 



^The association of 6i and 62 with derivatives of f{5m) is no longer accurate in the non-hnear regime, but one can 
still define 61 = [Pg{k) / Pm{k)]^^'^ and define 62 through equation jig), though 62 rnay now be dependent on triangle 
shape. 
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Fig. 12. — Reduced bispectrum Qeq{k) for equilateral configurations as a function of scale, for 
different HOD prescriptions. In panels (a)-(c), the dotted curve shows Qeq{k) for a model with 

power-law A^avg(M), Mmin = 2.8 x lO^^/i^^M©, a = 0.5, and Average P(iV|iVavg)- Other curves 
show Qcq{k) for models with different Mmin (panel a), different a (panel b), or different P(A^|A'^avg) 
(panel c). Panel (d) sgows models with different P(A^|A^avg) for a = 0.8. Thick solid curves show 
Qeq{k) for the unbiased mass distribution. 
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Fig. 13. — The reduced bispectrum Q{9) as a function of triangle shape, where ki = 2/c2 = 1/iMpc ^ 
(k = 27r/A). The four panels show the same HOD models as in Fig. O. 
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dependence on P(A^|A''avg) for higher a reflects the greater importance of the 1-halo and 2-halo 
terms when massive halos are highly occupied, since the 3-halo term is independent of P(A^|A''avg) 
but the other two are not. Figures [l^ and show a slight dependence of Q(^) on i-*(A^|A''avg) 
for a = 0.5 and a depression of Q(^) at all angles for the broader distributions when a = 0.8. 

We have found the influence of HOD parameters on the reduced bispectrum relatively difficult 
to interpret cleanly, for several reasons. First, HOD changes influence both the individual terms 

2-, and 3-halo) contributing to the bispectrum and the relative importance of those terms, 
and it is not always clear which effect is more significant. Second, the use of a Fourier space 
measure makes it more difficult to separate terms based on scale — the 1-halo contribution to the 
three-point function, for example, must vanish beyond the virial diameter of the largest halos, but 
it can still have an influence on the bispectrum down to fairly low k. Finally, while the influence of 
an HOD change on the bispectrum or power spectrum may be easy to guess, the effect on the ratio 
Q is often much harder to predict. The reduced bispectrum is a natural statistic for describing 
mass clustering, especially in the perturbative regime. However, the unreduced, conflguration 
space three-point function may prove more useful in unraveling the HOD, at least on scales of 
individual halos. 

We have also computed the hierarchical amplitudes 6*3 and 5*4 (related to the skewness and 
kurtosis, respectively), as a function of smoothing scale. These results are qualitatively similar to 
those for the reduced bispectrum shown in Figure so we do not show them here. 

Despite the complexities, three-point correlations clearly have information content that is not 
present in ^^(r), because of their dependence on higher moments of P{N\Nsn^ and their greater 
weighting of massive halos. The dependence of Q on Mmim «> and P(A^|A''avg) is signiflcantly 
different from the dependencies found for ^g(r) in § 3. We therefore expect three-point correlations 
to play an important role in constraining the HOD and breaking degeneracies between bias and 
cosmology. The suppression of Qg relative to Qm in HOD models with a < 1 (seen for all of the 
cases in Figure 12) may help to explain the low amplitude that Jing fc Borner (199^ ) find for 



the configuration-space three-point function of LCRS galaxies, which is reduced relative to the 
N-body prediction for the mass distribution by about a factor of two. 



4.3. Void Probability Function 

We now turn to a statistic which focuses on the lowest density regions, the voids in the galaxy 
distribution. The statistic we examine is the void probability function (VPF), and it is defined as 
the probability Po{R) that a randomly placed sphere of radius R contains no galaxies. The best 



observational study of the VPF to date is the analysis of the CfA2 redshift survey by Vogeley et 



al. (1994), though 2dFGRS and SDSS should yield results in the near future. Little &: Weinberg 



(1994 ) showed that the VPF is sensitive to the details of biasing, using environmental bias models. 
In the HOD context, the VPF should depend on the expected number of halos contained within a 
given sphere, as well as the expected number of galaxies per halo. The VPF does not distinguish 
between halos that contain one galaxy and those that contain multiple galaxies, since one galaxy 
alone prevents a region from being empty. Therefore, the VPF should be very sensitive to the halo 



mass regime Mmin ^ M < Mi, where most halos contain either one galaxy or no galaxies. Benson 
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(2001) has analytically studied count probabilities in the HOD framework, and Casas-Miranda 
let al. (2001 ) have studied halo counts in voids. These are steps towards an analytic description 
of the galaxy VPF under HOD bias. Here we provide an approximate model that is useful for 
understanding our numerical results. 



The number of galaxies expected in a volume V that has a density contrast 5 is 



roc J„ 

iN\V,S) = vl dM- 



iVavg(M), 



(17) 



where the halo mass function is conditional on 6 because halos are, on average, more abundant 
and more massive in high density regions than in low density regions. The amplitude of dn/dM is 
enhanced by the same factor (1 + 5) that affects the mass density. However, the shape of the mass 
function also depends on environment: in low density regions, gravitational growth of structure is 
suppressed, and the halo mass function shifts towards lower masses. This dependence of the halo 



mass function on density is the basis for the Mo &: White (1996 ) calculation of halo bias. For our 
present purposes, it is useful to write equation ( [l7| ) as 



iN\V,S) = Vil + S)fdM^ 



iVavg(M), 



(18) 



where 6l denotes the Lagrangian density contrast associated with the evolved non-linear density 
contrast 5. In equation (jl^), the dispersal of halos by the expansion of underdense regions is 
represented by (1 + 6), while the shift of the mass function is represented by 

For a Poisson distribution of mean (A^), the probability of having = is P{0) = exp(— (N)). 
The actual probability distribution of galaxy counts depends both on P(A^[A''avg) and on the 
probability distribution of halo counts; however, approximating the combined effect as Poisson 
provides a useful guide for understanding the VPF. In the Poisson approximation, we can write 
an expression for the VPF by integrating over S, 



P{N = 0\V) 



d5P{5\V){l + 5)exp 



.r f°° dn 

-V / dM 

Jo dM 



(19) 



where P{5\V) is the probability that a sphere of volume V has a mass density contrast 6. If we 
ignored the dependence of dn/dM on Sl, the exponential factor would simply be exp{—Vng), and 
the VPF would thus not depend on the HOD. For Vfig ^ 1, the exponential factor will not be too 
far from unity, so the VPF will depend mainly on P{S). However, at large R, where Vug ^ 1, the 
exponential suppression is usually large, and changes to A''avg(M) can potentially have very large 
effects on the VPF. 



Figure 14 shows VPFs for a variety of HOD models, as well as for the unbiased mass 



distribution. Since the VPF is very sensitive to the number density of particles in the distribution 
being studied, we have randomly sampled the mass distribution to have the same space density 
of particles as the galaxy distributions. All HOD models shown have a power-law A'avg(Af). 
Figure [I^ shows that changing Mmin from 1.4 to 2.8 (in units of W^^h~^MQ) has little effect but 
that raising it to 7 and again to 11.2 has a dramatic impact, especially at large R. This behavior 
has a natural explanation in the context of equation (p!9[). When Mmin is low, the integral inside 
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the exponential factor is never far from fig. However, when Mmin is large, the break in the halo 
mass function can shift below M^[^ in the lowest density regions, decreasing the integral and 
producing an enormous increase in the exponential factor. Figure ^b shows the effect on the 
VPF of varying a. HOD models with higher values of a have somewhat higher probability of 
containing voids, since they have lower A'avg at masses just above Mmin, but the dependence of 
the VPF on a is clearly weaker than the dependence on Mmin. Figure |l^ shows the effect of 
varying P(A''|A^avg)- Models with wider P(A^|A''avg) distributions have a slightly higher probability 
of containing voids because they sometimes produce empty halos in the halo mass regime where 
-^avg ~ 1 — 2. Finally, Figure |l4| d shows that the VPF is completely insensitive to the distribution 
of galaxies within halos. 

The VPF complements correlation statistics by responding sensitively to the low mass end of 
P{N\M) and especially to the cutoff mass Mmin- The VPF is a special case of count probabilities 
or of the closely related probability distribution function (PDF) of smoothed density fields. We 
have also measured these full PDFs, and they respond roughly as one would expect: models with 
more "bias" produce wider PDFs. However, we have been unable to draw simple associations 
between specific features of PDFs and specific properties of the HOD, so we have chosen not to 
present these results in detail. 



4.4. Pairwise Velocity Dispersion 

The statistics we have examined so far characterize the clustering of galaxies in real space. 
However, there is also important information stored in galaxy peculiar velocities, which are 
gravitationally induced by the underlying mass distribution. One of the statistics most widely 
used to characterize the galaxy velocity distribution is the pairwise radial velocity dispersion, 



^lir) = (|(vi - V2) • ri2p) - ((vi - V2) • ruf , (20) 



where vi and V2 are the velocities of two galaxies separated by ri2, and the average is over all 
galaxy pairs of separation r = |ri2|. This quantity can be estimated from the 2-dimensional 
correlation function ^(rp,7r) (e.g., pavis fc Peebles 1983| ), and its relatively low observed value 
provided one of the strongest early arguments against an 0^ = 1 cosmology with an unbiased 
galaxy distribution ( pavis et al. 1985| ). Recent observations give = 500 - 600km s"^ at 
r ~ l/i~^Mpc for typical optical samples, with a weak dependence on scale but a strong 
dependence on galaxy color or morphological type (Jing, Mo, & Borner 1998; Zehavi et al. 2001) 



Narayanan et al. (2000| ) found that cr^(r) is highly sensitive to the details of biasing, and here we 



wish to explore this sensitivity in the HOD framework of bias. 

On small scales, cr^(r) will be dominated by pairs of galaxies within the same halo. If we 
assume that the halo velocity distribution is isotropic and isothermal, then the mean pairwise 
radial dispersion of galaxy pairs in a halo of mass M is 20^0-2 (M), where a{M) is the halo's 
one-dimensional velocity dispersion and ay is the galaxy velocity bias parameter defined in § 2. 
The pairwise dispersion of the galaxy population is the average value of 2ala'^{M) weighted by the 
fraction of galaxy pairs of separation r contributed by halos of mass M. Reference to equation (|lT|) 
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Fig. 14. Void probability functions (VPF) for different HOD prescriptions. In each panel, the 
dotted curve shows the VPF for a model with power- law Afavg(M), Mmin = 2.8 X IO^^H-^Mq, 
a = 0.5, Average P(iV|iVavg), and A7 = 0. Other curves show VPFs for models with different 
-^min (panel a), different a (panel b), different P(A^|A''avg) (panel c), or different galaxy distributions 
within halos (panel d). Thick solid curves show the VPF for the unbiased mass distribution. 
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for S.ihi'f') shows that this average is 

-1 



dM 



dn {N{N-l)),, 



1 



dM 



F 



'I I 



(21) 



where we have assumed Cgi^) = Cihi'i') ^ 1- 



At large scales, each member of a galaxy pair resides in a separate halo, and the velocity 
of each galaxy is the center-of-mass velocity of its halo plus a random velocity of 1-d dispersion 
a^(T^(M). Since the random velocities are uncorrelated with the center-of-mass velocities and 
with each other, the mean pairwise dispersion of galaxy pairs in halos of mass Mi and M2 at 
separation r is 0"i2('^) + a^(T^(Mi) -|- a^(T^(M2), where o"^2('^) is the pairwise dispersion of halos. 
In the 2-halo regime, the pairwise dispersion of the galaxy population is the average of this 
quantity weighted by the number of pairs in halos of mass Mi, M2, which is proportional to 
(1 + 62(r))iVavg(Mi)iVavg(M2). Using equation (|) we have 



cr; 



v,2h 



(r) 



nlil+Cgir)) 



-1 



(1 + Cuir)) alcj\Mi) + a^a^(M2) + a^r 



dM.^N.. 
dMi 



avg 



(Ml) 



dM2^N.. 
dMo 



avg 



(M2) 



2a' 



where 



' ^ Jm„ 



dm 



is the mean halo dispersion weighted by galaxy number. 



Cj2(M) 



Nb 



dM^iV,,g {M)bh {M)a^ (M) 



(22) 



(23) 



(24) 



is the mean halo dispersion weighted by galaxy number and halo bias factor (and h is the galaxy 
bias factor of eq. B), and 
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(25) 



is the mean pairwise halo dispersion weighted by number of galaxy pairs. Since (Ti2{r) can 
in principle depend on both Mi and M2, it cannot generally be pulled out of the integral in 
equation (25), though at large scales it may be close to the value o"^iin('^) predicted for the mass 
by linear theory, since halos obtain their relative velocities from large scale density perturbations. 

For a more detailed discussion, we refer the reader to |Sheth et al. (2001 ), who give a much 
more comprehensive analytic treatment of the pairwise dispersion and mean streaming in the 
halo occupation framework and address a number of subtle issues (like non-isothermality) that 
we have glossed over here. For our purposes, the crucial features of equations (^l|) and (p2|) are 
the following. In the 1-halo regime, depends on and {N{N — 1))^^ , and different halo 
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masses contribute on different scales just as they do in Cihi^)- In the 2-halo regime, cT^(r) has 
one contribution that is proportional to and another that is independent of a^- At large 
separations, where Cg{r) is small, the first contribution becomes constant, so the change in cT^(r) 
is entirely due to (o"^2 (''))• 



Figure 15 shows the pairwise velocity dispersion for a variety of HOD models and for the 
unbiased mass distribution. Increasing Mmin (Fig- |l5|a) slightly increases cT^(r) because it forces 
galaxies into more massive halos with high velocity dispersions. However, increasing a (Fig. plsjb) 
has a much more dramatic effect because it redistributes galaxies preferentially towards the most 
massive halos. For all a, o"^(r) has a broad maximum at r l/i~^Mpc, where the highest mass 
halos contribute a significant number of pairs. Higher a models have higher C7^(r) even at large 
scales because of the continuing contribution of ((T^(M)) (eq. [^2|). Figure [l5| c confirms that 
(T^(r) is independent of P(A^|A'^avg) on large scales, where the 2-halo term dominates. However, on 
small scales P(A^|A'^a,vg) has an important impact. While making P(A^|A'^avg) broader boosts the 
correlation function (Fig. it drives a'^{r) down by allowing ih('') to be dominated by pairs 
in low mass halos. For example, the Average model never has pairs in halos with Ni^yg{AI) < 1, 
but the Poisson model does. Consequently, the small scale pairs in the Average model necessarily 
come from higher mass halos with larger velocity dispersions than in the Poisson model. 

Figure ^d shows the influence of on cT^(r), which is, as one might expect, very strong. 
The dotted line shows the extreme case of = 0, where (T^(r) comes only from the relative 
velocities of halos themselves. This contribution goes to zero in the 1-halo regime. It is usually 
sub-dominant even on large scales, but it dominates the change in C7^(r) because the internal 
dispersion contribution becomes constant. Increasing or reducing galaxy random velocities by a 
factor of two {oy = 2 or 0.5) changes cr^(r') by a large amount at all scales; note, however, that 
values of so far from unity are probably implausible on dynamical grounds. The dot-long-dash 
line in Figure ^ shows the impact of a central galaxy on a model with = 1; the pairwise 
dispersion drops slightly because the central galaxy moves with the halo mean velocity. 

In the absence of velocity bias, ay ~ 1, the pairwise dispersion is a good diagnostic of a, which 
controls the fraction of galaxies in high dispersion halos ( ping. Mo, Sz Borner (199^ ) reach a similar 
conclusion). It also constrains P(A^|A'^avg)) with a rather different signature from £,g{r). For the 
ACDM cosmology and HOD models considered here, matching the observed amplitude and weak 
scale-dependence of CTy{r) requires a ^ 0.5 to suppress the peak at r ~ l/i~^Mpc associated with 
massive halos, and sub-Poisson fluctuations at low M to prevent low mass halos from dragging 
down (T^(r) at small r. These are the same features of P(N\M) required to match ^g{r). The 
pairwise dispersion is also sensitive to and to the amplitude of the mass power spectrum, since 
these determine the scale of halo velocity dispersions and of relative halo velocities. 



4.5. Redshift Space Distortions 

The dispersion of galaxy peculiar velocities in collapsed regions smears structure along the 
line of sight in redshift space; it is this induced anisotropy that allows cr^(r) to be inferred from 
redshift maps without measuring galaxy peculiar velocities individually. On large scales, coherent 
flows into overdense regions and out of underdense regions amplify the contrast of structure along 
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Fig. 15. — Mean pairwise radial velocity dispersions as a function of the pair separation r for 
different HOD prescriptions. In each panel, the dotted curve shows for a model with power-law 
N^yg{M), Mmin = 2.8 X lO^^/i^^M©, a = 0.5, Average P{N\Na.yg), and a„ = 1. Other curves show 
ay for models with different Mmin (panel a), different a (panel b), different P(A^|A^avg) (panel c), or 
different velocity distributions within halos (panel d). Thick solid curves show for the unbiased 
mass distribution. The dot-long-dash curve in panel (d) shows a model with = 1 but the first, 
central galaxy in each halo forced to move at the halo's mean velocity. 
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the line of sight (Bargent &: Turner 1977 ; Kaiser 1987; for an excellent and comprehensive review 



see Hamilton 199 j ). In the linear regime, the linear bias approximation 5g = b6m, and the distant 
observer approximation, the redshift-space power spectrum P^{k) is enhanced by a factor 

over the real-space power spectrum P^{k), where /? = /(f^m) ~ ^m^/^ (Kaiser 1987). For 



models in which galaxy bias is a more complex function of local environment, equation (|2q ) still 
applies on large scales, with b corresponding to the large-scale asymptotic value of the rms bias 



factor (Berlind et al. 2001). However, non- linear effects, especially those associated with velocity 
dispersions in collapsed regions, influence P^{k) out to scales 27r/k ~ 50 — 100/i~^Mpc or more 
(Cole et al. 1994). Measurement of /? with redshift-space distortions therefore requires very large 



redshift surveys and accurate modeling of non-linear effects. The analysis of the 2dF galaxy 
redshift survey by Peacock et al. (2001) provides the best measurement to date. 



We compute power spectra of our mass and galaxy distributions by CIC-binning them onto a 
200^ grid and applying an FFT. For the redshift-space power spectrum, we take the line-of-sight 
direction to be a Cartesian axis of the simulation cube, implicitly assuming that the whole 



simulation volume is far enough away to satisfy the distant observer approximation. Figure 16 
shows P^{k)/P^{k) for a variety of HOD models and for the unbiased mass distribution. The 
general behavior of these curves is as expected: suppression of the redshift-space power spectrum 
on small scales (high k), where non-linear effects dominate, and enhancement at low k. However, 
it is clear that the largest scales probed by our 141.3/i"^Mpc simulation cube are not yet in the 
linear regime described by equation ([2^), since there is no clear plateau in P^{k)/P^{k) and the 
maximum ratio for the mass distribution is substantially lower than the value 1.37 implied for 
P = 0.30-6 ^ 0.48. 

In Figure [l6|a, raising Mmin depresses P^{k)/P^{k) at all k. At small scales, where 
P^{k)/P^{k) < 1, this behavior reflects the increased suppression of line-of-sight structure by 
velocity dispersions in high mass halos. At large scales, one expects higher M^[^ models to have 



lower P^ {k) / P^{k) because they have higher h (lower /?), but the effect seen in Figure 16a may 
still be dominated by the increased residual impact of small scale velocity dispersions. All four 
HOD models have a = 0.5, so relative to the mass distribution the high dispersion halos are 
down-weighted, and the large scale bias factor is slightly less than one. Thus, all four models have 
higher P^[k)/P^{k) than the mass distribution. Raising a reduces P^{k)/P^[k) by increasing 



velocity dispersions and the bias factor, as seen in Figure 16b. The a = 0.8 model has nearly the 
same pairwise dispersion as the mass (Fig. [T5|b) and a large scale bias factor slightly larger than 
one (Fig. ^); as one might therefore expect, its P^ {k) / P^{k) ratio matches the mass on small 
scales and is slightly lower on large scales. Figure ^6|c shows that the form of P{N\Ns_-^g) has no 
significant impact on the power spectrum distortion. 



Figure [iq d shows the impact of velocity bias within halos. Even though all four models have 
the same mean velocity field, their power spectrum ratios are different on all scales probed by 
the simulation. The comparison shows that the suppression of P^{k)/P^{k) on small scales is 
dominated by dispersions within virialized halos. However, even the a^, = model shows a drop in 
P^{k)/P^ {k) relative to the value at large scales, showing that the dispersion of the halo velocities 
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Fig. 16. — Ratios of the redshift to real space power spectra jp^ as a function of wavenumber 
k for different HOD prescriptions. In each panel, the dotted curve shows P^{k)/P^{k) for a model 
with power- law iVavg(M), Mmin = 2.8 X lO^^/j-^M©, a = 0.5, Average P(A^|A^avg), and a„ = 1. 
Other curves show P^{k)/P^'{k) for models with different Mmin (panel a), different o; (panel b), 
different P{N\Na;vg) (panel c), or different velocity distributions within halos (panel d). Thick solid 
curves show P^{k)/P^{k) for the unbiased mass distribution. 
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themselves also plays a role. Above all, Figure p^ d confirms previous indications that small scale 
galaxy velocity dispersions influence the redshift-space power spectrum out to remarkably large 
scales. 

Redshift-space distortions can also be quantified using the ratio of quadrupole and monopole 
moments of P^{k) ( pole et al. 1994 ). We have measured this ratio as well as P^{k)/P^{k) and 
find similar dependence on HOD parameters. 



Although our simulation is not large enough to show it, we expect on the basis of Berlind] 
et al. (200l|) that the power spectrum ratio would approach the value defined by /? = Q^/bp 



at sufficiently low k, where 6p is the asymptotic value of Pg{k)/Pm{k). Even with a large 
galaxy redshift survey, constraining f3 requires modeling the non-linear effects on redshift-space 
distortions, which is most often done using the Peacock & Dodds (1994) model of a linear theory 
distortion modulated by random velocities with a single value of the velocity dispersion. With 
galaxy redshift samples of the size expected from the 2dF and SDSS, the limitation of this model 
introduces a systematic error that is larger than the statistical uncertainty ( |Cole et al. 1995 ; 
Hatton fc Cole 1999| ). The HOD seems the natural framework in which to develop a more realistic 
and more accurate model of non-linear redshift-space distortions. The analytic studies of Seljak 
(200ID and [White (200i| ) provide the b asis for such a model, which should improve constraints on 



f^m and h from future analyses of larger data sets. 



5. Group Statistics 

The statistics that we have examined so far treat individual galaxies as the basic units of 
structure. A useful complementary strategy is to define groups of galaxies and measure their 
properties. We identify galaxy groups in real space using the same friends-of- friends algorithm 
that we used to find halos in the mass distribution. We use a linking length of 0.2 times the mean 
inter-galaxy separation, and we only retain groups containing four or more galaxies. Since our 
galaxy distributions have a space density of O.Ol/i'^Mpc^^, the linking length is approximately 
0.9/i~^Mpc. With this choice, most galaxies within the same halo are identified as belonging to 
the same group, but galaxies in separate halos are only occasionally linked. However, there are 
situations in which the galaxy populations of two halos are linked into a single group, and others 
where the population of one halo is split into two or more components, so the match between 
groups and halos is not exact. For a galaxy redshift survey there are additional complications 
including the need for different angular and redshift linking lengths to account for peculiar 
velocities, scaling linking lengths with distance to account for varying galaxy density in a 
flux-limited survey, and, in surveys that employ flber spectrographs, accounting for incompleteness 
imposed by minimum fiber separations. (For discussion of some of these issues, see, e.g., ooreJ 



Frenk, &: White 1993 and Nolthenius, Klypin, &: Primack 1994.) Here we consider the idealized 



case and leave the problem of design and performance of group-finding algorithms to studies that 
are tailored to specific data sets. 

The most basic group statistic is the multiplicity function, the abundance of galaxy groups 
as a function of group richness. This can be measured for groups identified in a redshift survey 
( Maia, da Costa, fc Latham 198^ ; Ramella, Geller, & Huchra 198£; Ramella et al. 1999), or from 
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the galaxy surface density in an imaging survey ( pott h Turner 1977 ; de Filippis et al. 1999D . We 



define the cumulative multiplicity function ngrp(> N) to be the number density of galaxy groups 
that contain or more galaxies. To the extent that groups and halos correspond one-to-one, 
the differential group multiplicity function is simply a convolution of the halo mass function with 
P[N\M)^ and the cumulative multiplicity function is 



>^) = EX (27) 



In the limit of a narrow P{N\M), the space density of groups of multiplicity is equal to the 
space density of halos with Ni^yg(M) = N. The direct connection implied by equation ([27|) makes 
the multiplicity function a promising tool for constraining P(N\M). 

Figure ^ plots ngrp(> N) for a variety of HOD models, starting in each case at the largest 
group in the box and moving up to higher densities as N gets smaller. The function becomes 



noisy at high N where the number of groups is small. Increasing Mmin (Fig. 17a) shifts ngj-p{> N) 
towards larger N, since the overall galaxy space density remains constant and the number of 
galaxies in high mass halos therefore increases. Increasing a has a much more drastic effect on 



the multiplicity function (Fig. 17b) because it boosts the multiplicity of higher mass halos by 



a larger factor, changing the slope of ngrp(> A'^). Broadening P(A^|A'avg) has a modest impact 



(Fig. 17c); while a given halo's multiplicity can scatter in either direction, there are more low 



mass halos to scatter to high N than vice versa, so the net effect is to boost the high-A^ tail of the 



multiplicity function. Finally, Figure |17| d shows that the group multiplicity function is completely 
insensitive to the spatial distribution of galaxies within halos, a reassuring indication that the 
friends-of-friends algorithm is consistently identifying the same groups regardless of their detailed 
internal structure. 

From the HOD perspective, the other promising aspect of groups is the possibility of 
measuring their masses dynamically and thus obtaining a direct association between A^ and M. 
For each group in our simulated galaxy population, we apply the projected virial mass estimator 
given by Heisler, Tremaine, &: Bahcall (19"85| ), 



where V^^i is the line-of-sight velocity of galaxy i relative to the group velocity centroid and R^^ij is 
the relative tangential separation between galaxies i and j. For this calculation, we treat one axis 
of the simulation box as the line-of-sight direction. (We have also tried the other mass estimators 
suggested by Heisler et al. 1985, with similar results.) 



Figure ^ plots group richness against estimated virial mass for a variety of HOD models. 
Figure |l8|a shows a model with power-law A''avg(M), a = 0.5, M^\^ = 2.8 x \^^'^h~^ Mq, and 
Average P(A^|A''avg)- Each small point represents an individual galaxy group, and large points 
show the mean richness in bins of Mvir. If the group-halo match and virial mass estimates were 
perfect, this figure would reproduce the top panel of Figure [l|, which illustrates the true P{N\M) 
used to create this galaxy population. However, random errors in Mvir make the scatter-plot 



in Figure 18 much broader. Perhaps more seriously, the values of Ng,yg{M) recovered from this 
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Fig. 17. — Cumulative group multiplicity functions for different HOD prescriptions. In each panel, 
the dotted curve shows ngi.p(> A^) for a model with power-law Na_yg(M), Mmin = 2.8 x 10^^/i~^Mq, 
a = 0.5, Average P(A^|A''avg), and A7 = 0. Other curves show ngrp(> N) for models with different 
Mmin (panel a), different a (panel b), different P{N\Na,yg) (panel c), or different galaxy distributions 
within halos (panel d). 
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diagram are biased with respect to the true relation, shown by the sohd hne. At M <^ lO^^/i~^M0, 
the mean multiphcity is overestimated because the group catalog excludes halos in this mass range 
that have < 4. At high masses, random errors in Mvir cause a slight underestimate of Nsfjg{M), 
since the number of low-A^ groups with overestimated Mvir exceeds the number of high-A^ groups 
with underestimated Mvir- Both of these biases are calculable, so data like those in Figure |l^ 
could be corrected to yield a more accurate estimate of A^avg(-M). However, such a correction 
would require some assumptions about P{N\M) (such as an extrapolation at the low mass end) 
and an accurate understanding of the error distribution of Mvir- 



The situation is considerably more promising for the a = 0.7 model shown in Figure [li 
This model has a much larger number of high multiplicity groups, which yield more accurate Myir 
estimates because they have more galaxy tracers. Furthermore, the dynamic range in A^ is large 
enough to allow an accurate estimate of the slope of A^avg(-M) at high M, though the amplitude of 
the relation is still slightly underestimated. Figure [l^ shows a model with a break in A^avg(-M) at 
Merit = I'S^^h-^MQ. The virial mass estimates reveal the existence of the break and recover the 
high-M slope of A^avg(M) reasonably well. However, the mass at which the break occurs and the 
slope at the low mass end of the relation are not clearly determined. A break at lower mass, like 
the one in the model illustrated in Figure ^, would pass undetected. 

Figure ^d shows an a = 0.5 model with a Poisson P(A^|A^a,vg)- Despite the scatter in the 
mass estimates themselves, the broader P{N\M) is clearly detectable in comparison to the Average 



model of Figure 18a. The larger scatter increases the bias in the recovered N.^^g{M) at low M. 
Figure I83 shows a model with A7 = +1, i.e., galaxies more extended than dark matter within 
halos, for which the recovered A^avg(-M) lies below the true relation at high M. Figure p!^ shows a 
model with a strong velocity bias, = 0.5, which drives estimated masses down relative to true 
values. For rich groups and clusters, mass estimates from galaxy dynamics can be tested against 
masses estimated from X-ray observations or gravitational lensing, and systematic errors as large 
as those in Figures p!8| e and [l8| f (which display fairly extreme models) can probably be ruled out 
on the basis of existing studies (see, e.g., |Wu et al. 1998|) . Such comparisons become harder for 
low multiplicity groups, though it may soon be possible to test for any overall bias by comparing 
mean dynamical masses to the group-mass correlation function derived from weak lensing. 

Figure suggests that it will be difficult to recover P{N\M) or N.^^g{M) simply by measuring 
dynamical masses as a function of multiplicity. Even if there are no internal biases (like ^ 1 
or A7 7^ 0) that systematically influence Mvir, the random errors in Mvir and absence of low 
multiplicity halos from the group catalog drive the recovered A^avg(M) away from the true value. 
However, viewed as another clustering statistic, the richness-virial mass relation can provide 
another valuable constraint on the HOD, especially if the form of A^a,vg(M) has been determined 
from other clustering measures, since it tests the overall mass scale and constrains the role of 
internal spatial and velocity biases. 

The poor performance of the mean richness as a function of virial mass in recovering A^avg(M) 
leads us to consider an alternative way of interpreting the richness virial mass relation. Figure |l9| 



shows this relation for the same set of HOD models as Figure 18, but large points now represent 
the mean virial mass at each multiplicity A^, and solid lines show the true Mavg(A^) relation. For 
a narrow P(A^|A^a,vg)) such as the Average distribution, Afavg(A^) and A^avg(M) contain identical 
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Fig. 18. — Richness vs. galaxy group virial mass estimates. Each panel represents a particular 
HOD model. The models are specified at the top of each panel, where M^nm is given in units of 
10^^h~^MQ, and Poi and Ave represent Poisson and Average P{N\Ng;vg) distributions, respectively. 

Galaxy groups were identified with a fricnds-of-friends algorithm in real space, and a virial estimator 
was used to estimate their mass. Each panel also shows the mean group richness in bins of Myir 
(large points) as well as the input Ni^^g(M) relation (line). 



-47- 



information, since A'avg(-^avg(-^)) = However, in the Poisson case these two functions differ, 



as can be seen in Figure |l9|d, where Ni^vg{M) is shown by the dashed hne. Figure [19| demonstrates 
that most of the biases that affected the recovery of Ng^^giM) do not affect the recovery of 
Ma.-vg{N). Only velocity bias seriously hampers the recovery of Mavg(A^) (Fig. p^), and, even in 
that case, only the amplitude of the relation is affected. Unfortunately, the Mavg(A^) relation 
cannot be inverted to recover Ni^^g{M) without assuming a form for P(A^|A'^avg)! so the improved 
performance of this statistic comes with the loss of model independence. 



6. Measuring the HOD Empirically 

The recurring theme of preceding sections is that different galaxy clustering statistics are 
sensitive to different aspects of the halo occupation distribution. The encouraging implication is 
that combinations of clustering measurements can pin down the properties of the HOD, hemming 
it in with complementary constraints. The 2dF and SDSS redshift surveys are ideally suited to 
this task because they can measure clustering with high precision over a wide dynamic range. 
The SDSS also provides detailed photometric information for all galaxies in the redshift survey, 
allowing studies of the HOD for a multitude of galaxy classes, and weak lensing measurements 
from the SDSS images can determine the galaxy-mass correlation function for each of these classes. 

In this section, we outline a strategy for determining the HOD from observations, and 
although this outline is far from complete, we believe that it offers a promising route forward. 
Since the possible parameter space of HOD models is extremely broad, a practical strategy 
requires some procedure for getting a good first guess at P{N\AI), which can then be refined 
and tested using a battery of clustering measures. Our procedure assumes that the underlying 
cosmological model is known in advance based on independent observations, though the model 
can be further tested by seeing whether it reproduces all aspects of observed galaxy clustering for 
some HOD. The ability of large scale structure to constrain cosmology depends in part on whether 
two different cosmological models with two different halo occupation distributions can produce the 
same galaxy clustering; we consider this question briefly at the end of the section. 

We take as our starting point the tight connection between P{N\M) and the group multiplicity 
function implied by equation (p7|). If we assume that P(A^|A'^avg) is narrow, then group multiplicity 
is a monotonic function of halo mass, and determining Ns;yg(M) from the multiplicity function is 
simply a matter of matching space densities. Given the measured ngi.p(> N) and the cumulative 
halo mass function nh{> M) computed from the cosmological model, we choose Ni^-vg{M) so that 
%rp[> -^avg(-^)] = nh{> M). More generally, one can assume a form of P(A^|A''avg) and infer 
-^avg(Af) by reversing the convolution in equation (p7|). The relatively small impact of P(A'"|A^avg) 



on the multiplicity function compared to changes in A'^avg(Af) (see Fig. 17) shows that results will 
not be overly sensitive to the assumed width of the distribution. 



Figure 20 illustrates the performance of the monotonic matching method. The top panel plots 
the cumulative halo mass function derived from the GIF N-body simulation, with halos identified 
by the FoF algorithm as usual. The middle panel shows cumulative group multiplicity functions 
for three HOD models, with groups identified by the same FoF algorithm as described in § 5. The 
bottom panel shows the A''avg(M) relations derived by matching nh{> M) to ng^p{> N) and the 
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Fig. 19. — Same as Fig. 18, except that large points represent the mean estimated virial mass at 
fixed group richness N, sohd hnes show the true Ma_yg(N), and the dashed hne in panel (d) shows 

iVavg(M). 
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true Ni^^g{M) used as inputs to the HOD models. In all three cases, the method is impressively 
successful at recovering N^vgiM) for M ^ W^^'^h~^ Mq, obtaining the correct power-law slopes 
for the a = 0.6 and a = 1.0 models and revealing the break in the broken power-law model. At 
low multiplicities, the discreteness of N causes a step-like behavior in Ni^^g{M), but this could be 
easily repaired. The method works well for the a = 1.0 Poisson model even though we have not 
applied any deconvolution correction. 



The matching method illustrated in Figure 2C is a promising way of determining N^vgiM) at 



high M, where halos contain large numbers of galaxies. For application to redshift survey data, the 
main challenge is to develop a good group-finding algorithm and understand any biases associated 
with it. When studying the HOD of galaxy classes, one would start with the group catalog and 
associated halo masses derived from the full, high density galaxy catalog, then measure the content 
of individual classes. 

Group matching provides a partial guess at the HOD, which must be extrapolated to 
low masses, refined, and tested using other statistics. One important constraint on the 
extrapolation comes from the mean galaxy number density, fig = dM-^Ns;^g{M). A 
similar constraint with a different weighting of halos comes from the large scale bias factor, 
h = dM ■j^bh{M)Nscvg{M), which can be inferred from the measured galaxy power 

spectrum given the assumed cosmological model. The most valuable constraints on halo 
occupation at low masses comes from the 1-halo regime of ig{r), since different separations probe 
different halo mass scales. 

We can illustrate the constraining power of ^g{r) using an approximation suggested by 
Figure which shows that the cumulative pair separation distribution F(r/2i?vir) is roughly a 
linear ramp out to a maximum separation Xm ~ 0.6 of the virial diameter, and is close to unity at 
larger separations. Under this approximation, the 1-halo term ■^ih('") only receives contributions 
from halos with virial diameters 2i?vir > , or, using equation (|l^), with virial masses 

M,(r) = '^^{j-)\ (29) 

We can therefore rewrite equation ( [ll| ) by substituting the constant value F'{x) = x^ and 
changing the lower limit of the integral to Mi{r), which is now the only part of the integral with 
an r dependence. Differentiating both sides of the equation with respect to r and rearranging 
terms yields an expression relating the second factorial moment at Mi to the 1-halo correlation 
function and its logarithmic derivative at r: 



{N{N-1)),,^ 




dlnM 



M, 



3 , (iKih(r) \ 
^ ( 2 + ) 6h(r). (30) 



Figure 21 illustrates the performance of equation (|30|) . Despite the rather crude, linear ramp 
approximation for F[r/2R^i^), the estimates of {N{N — 1))^ derived from Cih(?") recover the true 
values shown by the solid lines for both of these HOD models, which have Mmin = 2.8 x lO^^/i^^M©, 
a = 0.5, and Poisson and Average P(A^|A'avg), respectively. In practice, one could only apply 
this method to halos with virial diameters 2i?vir ^ O.Sx^^/i^^Mpc ~ 0.8/i^^Mpc (mass 
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Fig. 20. — Recovery of A''avg(M) from the group multiplicity function. Panel (a) shows the halo 
mass function nh(M) derived from the GIF N-body simulation, where halos were identified with 
a FoF algorithm. Panel (b) shows group multiplicity functions ngip(> A'') for three HOD models, 
where galaxy groups were identified with the same FoF algorithm. Panel (c) shows the resulting 
^avg{M) relations from matching nh(M) to ngj.p(> N) (points). Also shown, for comparison, are 
the true Na.yg(M) functions that were used as inputs for the HOD models (lines). 
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Af ^ 5 X for = 0.3), since ^g(r) departs from ^ih(r) at r ^ 0.5/i-^Mpc (see Fig. 0). 

However, one could extend the range somewhat by calculating the 2-halo contribution to (,g{r) and 
subtracting it to obtain Cihi"!"); alternatively, a numerical method that uses equation (11) directly 
instead of the linear ramp approximation would not be difficult to implement. The important 



lesson of equation ( |30D and Figure 21 is that ^ih(?') provides fairly well localized information about 



{N{N — 1))^,/, in just the mass range where the group matching method breaks down. 

The small scale correlation function depends on {N{N — l))j^^ rather than A'avg(Af), and 
it has an additional dependence on the assumed galaxy profile (though Fig. |6| shows that this 
dependence is not strong). However, it is clear from the discussion above that matching ngi-p 
Ug, and (,g{r) already places strong constraints on P(N\M) for a given cosmology. A trial HOD 
inferred by the techniques outlined above can be tested and refined using other clustering statistics. 
The VPF gives strong guidance on the low mass end of P{N\M), especially the value of Mmin. 
Higher order correlations probe the higher moments of P^NlNgj^g). Additional strong tests come 
from the mean virial mass as a function of multiplicity and from the pairwise velocity dispersion, 
which is sensitive to the mean occupation at high mass and to -P(iV|iVavg) at low mass. Success in 
matching spatial clustering but failure with these dynamical measures could be a signature of an 
incorrect assumption about the underlying cosmology (e.g., the wrong value of Qm), or it could be 
an indication of velocity bias (a^ 7^ 1). The scale dependence of the pairwise dispersion can help 
discriminate between these possibilities, since cr^(r) combines terms that depend on with others 
that do not (see eq. |2^). Lensing measurements of the galaxy-mass correlation function (or the 
analogous group-mass correlation function) provide stronger discrimination by probing halo mass 
scales independently of galaxy peculiar velocities. The 1-halo regime of Cgmif) also constrains 
in a range where (.gii") constrains {N(^N — 1))^^^. 

Suppose we find an HOD that reproduces all of these properties of observed galaxy clustering. 
Can we take this success as evidence that the assumed cosmological model is itself correct, 
or could there be some other combination of cosmology and HOD that would yield the same 
results? In other words, do cosmology and HOD bias have degenerate effects on galaxy clustering, 
or can sufficiently precise measurements constrain both independently? We reserve a detailed 
investigation of this question for future work, but qualitative arguments suggest grounds for 
optimism. 

Let "model 1" denote our supposed successful combination of cosmological model and HOD 
prescription. If we raise fi^ but keep the linear theory matter power spectrum P\i^{k) the same, 
then the dark matter halo population will be essentially unchanged except that all masses will 
grow by the factor ^lm,2/^m,i, including the characteristic mass M* of the halo mass function. If 
we maintain the same HOD as a function of M/M*, then the spatial clustering of galaxies will 
stay the same, but dynamical measures like the pairwise velocity dispersion or virial mass vs. 
richness relation will not. Changes to dynamically sensitive statistics could be partially repaired 
by introducing velocity bias, with ay^2 < Civ,i, but this suppression will not disguise changes in 
cr^(r) or P''^{k) at large scales or changes in ^gm{r) at any scale. 

For a less easily detectable change, we can reduce the amplitude of Pim{k) at the same time 
we increase ftm, maintaining a "cluster normalization" condition as,2^^2 ~ ^8,i^m^i ( [White, 



Efstathiou, &: Frenk 1993). This change preserves (approximately) the velocity scale of rich clusters 
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Fig. 21. — Recovery of {N{N — from ^ih{r) using eq. (^0[). The solid and open circles show 
the recovered {N{N — 1))^^^ for Poisson and Average -P(A^|A''avg) HOD models, respectively. Also 
shown for comparison, are the "true" {N{N — relations that were used as inputs for the HOD 
models (lines). 
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and large scale flows, but it reduces the length scale of non-linearity in the matter distribution, by 
roughly the factor {i^m,2/^m,i)^^^ required to keep constant in physical units. If we maintain 
the same relative galaxy populations as a function of M/M^,, then all characteristic scales in the 
galaxy clustering drop by the same factor, including the galaxy correlation length rg. To return 
ro and the space density fig to their original values, we would need to raise Mmin/-^*, or change 
a, or make other alterations to P{N\M), and these changes would generally have different impact 
on different clustering statistics. Furthermore, the cluster normalization condition preserves the 
amplitude of the halo mass function on cluster scales but changes the overall shape of dn/dM in 
physical units, with a steeper fall-off at high masses for models with lower erg (higher Vtm)- This 
change could be partially repaired by altering the shape of P\\n{k)-, but that would alter the halo 
spatial clustering. We tentatively conclude that substantial changes to the value of or the 
normalization or shape of P\\a{k) cannot be masked by changes in the HOD. (See [Zheng et al] 



1 2002 1 for further discussion of these issues.) 



Given the increasingly tight constraints on cosmological parameters that come from CMB 
anisotropics, the Lyman-a forest, gravitational lensing. Type la supernovae, cluster evolution, and 
other data, a procedure that recovers the HOD for an assumed cosmology would be very valuable 
in itself. Successful application of this procedure to surveys like the 2dFGRS and SDSS would tell 
us essentially everything that local {z = 0) galaxy clustering has to say about galaxy formation, 
assuming that the adopted cosmological model is correct. Determinations of the HOD for different 
galaxy classes would provide detailed targets for analytic or numerical models of galaxy formation 
and insight into the physics that controls galaxy luminosities, colors, and morphologies. 

If, as we have speculated, the effects of cosmology and bias are not degenerate, then the 
HOD determination procedure becomes something still more powerful: a systematic way of testing 
cosmological models against galaxy clustering data given only some very general assumptions about 
galaxy formation (basically the validity of the HOD framework itself). Incorrect cosmological 
models would fail to reproduce the observations for any choice of HOD. Of course, galaxy 
clustering provides a testing ground for cosmology even with less detailed assumptions about 
galaxy formation, since the effects of "local" bias on the power spectrum, higher order moments, 
redshift-space distortions, and galaxy density and velocity fields become relatively simple on 



sufficiently large scales (see, e.g.. 


Coles 19931; Fry & Gaztahaga 1993|; 


Mann, Peacock, & Heavens 


1998 


; Bcherrer & Weinberg 1998 


; Narayanan et al. 2000|; 


Berlind et al. 2001). However, by 



incorporating HOD determination into analyses of galaxy clustering, one can use information from 
the highly nonlinear scales of individual halos to calculate effective "bias factors" for a variety of 
clustering statistics in the intermediate, moderately nonlinear regime, where the impact of bias 
may not be simple. This approach sharpens the sensitivity of clustering studies to small departures 
from the "minimal" cosmological model, which could be produced by light but not massless 
neutrinos, by a non-standard relativistic background, by non-scale invariant inflation, or even by 
a small contribution from isocurvature or non-Gaussian fluctuations. Analysis of separate galaxy 
classes with different HODs would provide redundant checks of the results, since the properties of 
the underlying mass clustering should be independent of the tracer population used to infer them. 

The monotonic matching method proposed above is similar to the technique suggested by 
Peacock &: Smith (2000| ), except that they use the "galaxy systems luminosity function" in place 



of the multiplicity function and extend the assumption of a monotonic relation between halo 
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mass and luminosity into the single-galaxy regime. Applying the technique to the data of Moore 
et al. (1993) and the ACDM halo mass function, Peacock &: Smith (2000| ) recover an A''avg(Af) 
relation that is shallow at low mass and steepens towards higher masses, qualitatively similar 
to the predictions of semi-analytic models and hydrodynamic simulations ( [Benson et al. 2000 ; 
Berlind et al., in preparation) and to the broken power-law model whose correlation function is 
illustrated by the solid curve in Figure ^. As the first determination of the HOD from empirical 
data, this recovery of a relation that resembles theoretical predictions and yields approximately 
the observed galaxy power spectrum ( Peacock &: Smith 2000| , figure 8) is very encouraging. 
Using similar methodology, Marinoni &: Hudson (2001 ) have recently obtained a similar Ni^^g^M) 
relation for a new sample of optically selected galaxies. With a different, model fitting approach, 
Jing, Mo, &: Borner (1998D , |Jing fc Borner (199^) , and |Jing, Borner fc Suto (20021) conclude that 
suppression of the fraction of galaxies in high mass halos can help explain the observed form and 
amplitude of the galaxy two-point and three-point correlation functions and pairwise velocity 
dispersion. Determinations of the HOD from the 2dF and SDSS redshift surveys will yield a far 
more detailed picture of the relation between galaxies and dark matter and among different galaxy 
types themselves. The improved constraints on dark matter clustering have the potential to reveal 
effects that are quantitatively subtle but have profound physical implications. 



7. Summary 

In the HOD framework, the bias of a population of galaxies is fully defined by the conditional 
probability P{N\M) that a halo of virial mass M contains galaxies, together with prescriptions 
that specify the relative spatial and velocity distributions of galaxies and dark matter within 
virialized dark matter halos. This framework provides an informative basis for connecting 
observations of galaxy clustering to the physics of galaxy formation, and it opens a route to the 
empirical determination of bias from galaxy redshift surveys. We have investigated the sensitivity 
of galaxy clustering statistics to features of the HOD, focusing on models with mean occupation 
A''avg(^) oc for M > Mmin, and found the following results: 

1) The amplitude of the galaxy correlation function on large scales is determined by the 
relative numbers of galaxies in high mass and low mass halos, and is thus sensitive to a and, to a 
lesser extent, Mmin. On small scales, Cgi^) receives the largest contribution from galaxy pairs in 
halos of virial radius i?vir ~ t, so different separations probe the occupation at different masses. 
The amplitude and shape of (,g{r) are highly sensitive to a, Mmin, and P(A^|Aavg), though each 
feature affects (,g{r) in a different way. The spatial distribution of galaxies within halos has a 
relatively modest impact that is confined to small spatial scales. 

2) In HOD models, a power-law correlation function emerges from a balance of several 
competing effects, and it therefore requires rather specific combinations of parameters. The 
sensitivity of ^g(r) to HOD parameters implies that its observed power-law form is a strong 
constraint on the physics of galaxy formation, and that the success of semi-analytic models and 
hydrodynamic simulations in reproducing this form is entirely non-trivial. 

3) The influence of HOD parameters on the galaxy-mass cross-correlation function is similar 
to the effect on Cg(r), except that (,gm{r) is independent of i-'(A^|A''avg) even on small scales. The 
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direct dependence of galaxy-galaxy lensing measurements on the mass of halos makes them an 
especially valuable test of the background cosmology. 

4) The bispectrum contains complementary information to two-point statistics due to its 
greater weighting of high mass halos and its dependence on the third moment of P(A^jA'^avg)- The 
dependence of the reduced bispectrum Q{k) on a and M^ain is fairly complex, and different from 
the dependencies found for ^g(r). 

5) The void probability function (VPF) is sensitive to the low M regime of P{N\M), where 
halos have a significant probability of containing no galaxies. In particular, the VPF is most 
sensitive to the cut-off halo mass Mmin. 

6) The pairwise velocity dispersion is particularly sensitive to a, which controls the fraction of 
galaxies in high dispersion halos, and to the presence of velocity bias in halos. In addition, <T^(r) 
is sensitive to $7„t and the amplitude of the matter power spectrum. 

7) The ratio of the redshift to real-space power spectra is sensitive to Na_yg(M) and velocity 
bias, as well as Jl^. Velocity dispersions within halos influence the redshift-space power spectrum 
out to large scales. 

8) The group multiplicity function is sensitive to Mmin and a in different ways: increases in 
Mmin cause a roughly horizontal shift of ngrp(> N) towards greater N, whereas increases in a 
cause a drastic change in the slope of ngrp(> N). 

9) The richness - virial mass relation is a direct probe of P{N\M). However, it is subject to 
biases that make it difficult to extract N^^-vg^M). The mean halo mass per multiplicity, Mavg(A^), 
can be extracted much more robustly, though even this measure can be affected by velocity bias. 

The different sensitivities of these clustering statistics to different features of the HOD suggest 
that the HOD can be empirically determined from observations of galaxy clustering. We have 
outlined a strategy for measuring the HOD from a galaxy redshift survey, which involves obtaining 
a first guess at P{N\M) and then refining that guess with additional clustering statistics. In its 
present form, our strategy assumes that the cosmological model is known a priori. However, the 
effects of changing the HOD should generally be different from the effects of changing cosmological 
parameters, so an incorrect cosmological model may be unable to reproduce the full suite of 
galaxy clustering statistics for any choice of HOD. At the very least, determinations of the HOD 
for different classes of galaxies from the 2dF and SDSS redshift surveys will provide considerable 
insight into the physics of galaxy formation and the origin of galaxy properties. If the degeneracy 
between cosmology and bias can indeed be broken, then determinations of the HOD will also 
sharpen tests of the standard cosmological model, perhaps revealing the imprint of new physics on 
the large scale structure of the universe. 
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